Solving the 1-D Wave Equation with RNPL

The goal of this tutorial is quickly guide you through the use of a pre-coded RNPL application that solves a simple time-dependent PDE using finite difference techniques. You can use the code from this example as a template for your work in Project 1.

The problem that we consider is the one-dimensional wave equation:

phi(x,t)_tt = phi_xx

where we again adopt the subscript notation for partial differentiation, and we have set the speed of wave propagation to unity. As discussed in class, the general solution of this equation represents a superposition of left moving (l) and right moving (r) waves (radiation):

phi(x,t) = l(x + t) + r(x - t)

We will solve the equation on the domain {x = [0..1] : t = [0..T]} with outgoing radiation or Sommerfeld conditions specified at x = 0 and x = 1. That is, at x = 0 and x = 1, we demand that there be no radiation coming from the left or right respectively. Thus we have

phi_t = phi_x at x = 0

and

phi_t = -phi_x at x = 1

As also discussed in class, we recast the wave equation in first order form (first order in time, first order in space), by introducing auxiliary variables, pp and pi, which are the spatial and temporal derivatives, respectively, of phi:

pp(x,t) = phi_x

pi(x,t) = phi_t

The wave equation then becomes the following pair of first order equations

pp_t = pi_x

pi_t = pp_x

and the boundary conditions are

pp_t = pp_x at x = 0

pp_t = -pp_x at x = 1

pi_t = pi_x at x = 0

pi_t = -pi_x at x = 1

The RNPL code w1dcn_rnpl, solves the above system of equations using a second-order Crank-Nicholson approximation: second-order centred approximations to the spatial derivatives are used for the interior equations, while second-order forward and backward formulas are used for the spatial derivatives needed in the implementation of the left and right boundary conditions respectively.

We will work through the RNPL code in some detail in the tutorials, and you can also refer to the online documentation for the language for additional information.

Building the RNPL application

A tar-ed and gzip-ed distribution for the RNPL solution of the above wave equation is available for download HERE. I suggest that you create a new directory pi-nr, move the distribution file to that directory then execute

% tar xzf w1dcn.tar.gz

Change to the w1dcn directory and generate a listing of the files therein. You should see the following

% cd w1dcn
% ls
CVS  id0  id1  id2  id3  id4  Makefile  w1dcn_rnpl

Build the application by typing make; the make command should generate the following output.

% make
rnpl -l allf  w1dcn_rnpl
rnpl_fix_f77 updates.f initializers.f residuals.f
gfortran -fno-second-underscore -O6 -c w1dcn.f
gfortran -fno-second-underscore -O6 -c updates.f
gfortran -fno-second-underscore -O6 -c residuals.f
gfortran -fno-second-underscore -O6 -fno-second-underscore -L/usr/local/lib w1dcn.o updates.o residuals.o -lrnpl -lxvs -o w1dcn
gfortran -fno-second-underscore -O6 -c w1dcn_init.f
gfortran -fno-second-underscore -O6 -c initializers.f
gfortran -fno-second-underscore -O6 -fno-second-underscore -L/usr/local/lib w1dcn_init.o updates.o residuals.o initializers.o -lrnpl -lxvs -o w1dcn_init

Note that the rnpl compilation generates a number of Fortran source files (with .f extensions) and include files (source code fragments with .inc extensions)

% ls *.f *.inc
gfuni0.inc   initializers.f  residuals.f    updates.f  w1dcn_init.f
globals.inc  other_glbs.inc  sys_param.inc  w1dcn.

as well as two executables, w1dcn, which performs the evolution of the FD system, and w1dcn_init, which generates the intial data for the evolution:

% ls w1dcn w1dcn_init
w1dcn*  w1dcn_init*

If you are so inclined, you may wish to quickly inspect the contents of initializers.f, residuals.f and updates.f.

Running the RNPL application

To run the RNPL application, simply invoke w1dcn with a parameter file supplied as the single command line argument. For example, to run with the level-0 (base level discretization) parameter file, id0, type

% w1dcn id0
 Can't open in0.sdf
 Calling initial data generator.
 WARNING: using default for parameter epsiterid.
 WARNING: using default for parameter maxstep.
 WARNING: using default for parameter maxstepid.
 WARNING: using default for parameter s_step.
 WARNING: using default for parameter ser.
 WARNING: using default for parameter start_t.
 WARNING: using default for parameter trace.
 WARNING: using default for parameter epsiterid.
 WARNING: using default for parameter maxstep.
 WARNING: using default for parameter maxstepid.
 WARNING: using default for parameter s_step.
 WARNING: using default for parameter ser.
 WARNING: using default for parameter start_t.
 WARNING: using default for parameter trace.
Starting evolution.  step:            0  at t=   0.0000000000000000
 step:            1  t=  1.25000000000000007E-002  steps=           7
 step:            2  t=  2.50000000000000014E-002  steps=           7
 step:            3  t=  3.75000000000000056E-002  steps=           6
 step:            4  t=  5.00000000000000028E-002  steps=           7
 step:            5  t=  6.25000000000000000E-002  steps=           7
 step:            6  t=  7.49999999999999972E-002  steps=           7
 step:            7  t=  8.74999999999999944E-002  steps=           7
                            .
                            .
                            .
 step:          123  t=   1.5374999999999965       steps=          11
 step:          124  t=   1.5499999999999965       steps=          11
 step:          125  t=   1.5624999999999964       steps=          11
 step:          126  t=   1.5749999999999964       steps=          11
 step:          127  t=   1.5874999999999964       steps=          11
 step:          128  t=   1.5999999999999963       steps=          11

The various diagnostic messages written to the terminal require some explanation. First,

Can't open in0.sdf
Calling initial data generator.

indicates that the intial state file, in0.sdf, that is specified in the parameter file id0 does not exist in the current directory, and thus the initial data generator, w1dcn_init, is automatically invoked to prepare it. w1dcn_init is passed the same parameter file used in the invocation of w1dcn, and then goes about its job of generating the initial state file. Once w1dcn_init has completed this task, w1dcn resumes execution.

During the initialization phases of both w1dcn_init and w1dcn, warning messages are issued if default parameter values are being used for any of the generic RNPL parameters (generic in the sense that they are defined for all RNPL applications). Thus, in the diagnostics, we see the sequence

WARNING: using default for parameter epsiterid.
WARNING: using default for parameter maxstep.
WARNING: using default for parameter maxstepid.
WARNING: using default for parameter s_step.
WARNING: using default for parameter ser.
WARNING: using default for parameter start_t.
WARNING: using default for parameter trace.

appearing twice, first from the execution of w1dcn_init, then from the execution of w1dcn. These particular messages can be safely ignored.

The final sequence of messages:

Starting evolution.  step:            0  at t=   0.0000000000000000
 step:            1  t=  1.25000000000000007E-002  steps=           7
 step:            2  t=  2.50000000000000014E-002  steps=           7
 step:            3  t=  3.75000000000000056E-002  steps=           6
 step:            4  t=  5.00000000000000028E-002  steps=           7
 step:            5  t=  6.25000000000000000E-002  steps=           7
 step:            6  t=  7.49999999999999972E-002  steps=           7
 step:            7  t=  8.74999999999999944E-002  steps=           7
                            .
                            .
                            .
 step:          123  t=   1.5374999999999965       steps=          11
 step:          124  t=   1.5499999999999965       steps=          11
 step:          125  t=   1.5624999999999964       steps=          11
 step:          126  t=   1.5749999999999964       steps=          11
 step:          127  t=   1.5874999999999964       steps=          11
 step:          128  t=   1.5999999999999963       steps=          11

represents a trace of the time-stepping alogorithm of the application. The time steps at which diagnostic output, such as

 step:            1  t=  1.25000000000000007E-002  steps=           7

is produced, coincide with the those at which grid function output to .sdf files is generated. Both are controlled by the output parameter, which in id0 is set via

output := 0-*

indicating that output (at level 0) is to occur at every time step. In addition to the step number and physical time, the steps = ... diagnostic shows the number of sub-iterations that were required to achieve convergence of the difference equations at that particular time step. The convergence criterion is controlled by the parameter epsiter, which is set in id0 via

epsiter := 1.0e-5

Grid function output (.sdf files)

Grid function output from an RNPL application is stored in special, binary (i.e. non-human-readable), machine-independent files, known as .sdf files (for "Scientific Data Format"). By convention, such files have an .sdf extension, and RNPL will always append the string "_[level]" to the base part of the .sdf filenames that it generates. Thus, for example, our execution of w1dcn id0 produces two files with names ending in "..._0.sdf" (i.e. level-0 output files)

% ls *_0.sdf
pi_0.sdf  pp_0.sdf

one for each of the grid functions for which output has been enabled in the RNPL source code via the {out_gf = 1} directive:

float pp on g1 at 0,1 {out_gf = 1}
float pi on g1 at 0,1 {out_gf = 1}

The sdfinfo command can be used to get summary information about the contents of an .sdf file:

but more typically, and as discussed briefly below, the data stored in an .sdf file is sent to a visualization utility, such as xvs or DV for inspection and/or analysis.

Visualizing and analyzing 1D RNPL output using xvs

Refer to the documentation HERE for more complete information on the use of xvs and related utilties.

Start xvs from the command line as follows:

% xvs &
xvs: Starting service on port 5001

xvs: For help/documentation see

http://laplace.physics.ubc.ca/Doc/xvs/

(your invocation of the command will use a different port number)

Once xvs is running, you can transmit data to it from 1-D .sdf files using the sdftoxvs command. For example, if you invoke

% sdftoxvs pp_0.sdf

or, equivalently

% sdftoxvs pp_0

the xvs interface should then look like this:

You can then use the utility to animate the data, step through it one frame at a time, perform manipulations such as taking numerical derivatives, and much more---again see the available on-line documentation for more information.

Convergence testing the RNPL application

RNPL was designed to make convergence testing easy. The w1dcn distribution includes 5 initial data files, id0, id1, id2, id3 and id4. As can be readily verified using the diff command, the only differences in the parameter settings in these files are for level, in_file and out_file:

% diff id0 id1
8c8
< level := 0
---
> level := 1
13,14c13,14
< in_file  := "in0.sdf"
< out_file := "out0.sdf"
---
> in_file  := "in1.sdf"
> out_file := "out1.sdf"

In particular, all 5 id files define the same solution domain and initial data parameters. By running w1dcn with each of the parameter files, we generate a sequence of approximate solutions to the same continuum problem, but with grid resolutions of h (level 0), h/2, h/4, h/8 and h/16. For higher-level runs (i.e. with level other than 0), the RNPL code automatically adjusts the frequency of output so that identical physical output times result from each run. For example, performing a level-1 run, we see

% w1dcn id1
 Can't open in1.sdf
 Calling initial data generator.
 WARNING: using default for parameter epsiterid.
 WARNING: using default for parameter maxstep.
 WARNING: using default for parameter maxstepid.
 WARNING: using default for parameter s_step.
 WARNING: using default for parameter ser.
 WARNING: using default for parameter start_t.
 WARNING: using default for parameter trace.
 WARNING: using default for parameter epsiterid.
 WARNING: using default for parameter maxstep.
 WARNING: using default for parameter maxstepid.
 WARNING: using default for parameter s_step.
 WARNING: using default for parameter ser.
 WARNING: using default for parameter start_t.
 WARNING: using default for parameter trace.
 Starting evolution.  step:            0  at t=   0.0000000000000000
 step:            2  t=  1.25000000000000007E-002  steps=           5
 step:            4  t=  2.50000000000000014E-002  steps=           5
 step:            6  t=  3.74999999999999986E-002  steps=           5
 step:            8  t=  4.99999999999999958E-002  steps=           5
 step:           10  t=  6.24999999999999931E-002  steps=           5
 step:           12  t=  7.49999999999999972E-002  steps=           5
             .
             .
             .
 step:          246  t=   1.5375000000000050       steps=          10
 step:          248  t=   1.5500000000000052       steps=          11
 step:          250  t=   1.5625000000000053       steps=          11
 step:          252  t=   1.5750000000000055       steps=          11
 step:          254  t=   1.5875000000000057       steps=          11
 step:          256  t=   1.6000000000000059       steps=          10

Note that output (both diagnostic and .sdf) is now generated every two steps, but that the physical output times are commensurate with those from the level-0 run.

We can quickly perform a series of runs at levels 0 through 4 (output supressed)

% w1dcn id0
% w1dcn id1
% w1dcn id2
% w1dcn id3
% w1dcn id4

then send all of the data for pp to xvs

% xvs ka  Closes pre-existing xvs windows
% sdftoxvs pp_*.sdf

The xvs display should now look like this:

We can then have xvs check the convergence of pp by selecting O(h^2) convergence test all from the applyn xy sub-item of the pull down menu (bring up the menu with the right mouse button). The resulting computation will generate another xvs window, leaving the display looking as follows:

Bring the window labelled pp_4_ to the foreground by positioning the cursor in it and hitting the space bar. Type S, then use the right and left mouse buttons to step forward and backward, respectively, through the results of the convergence test. For example, stepping to the 28th frame, you should see

Each frame (time step) of the data shown in the pp_4_ window displays 4 distinct quantities, namely

       pp_1 - pp_0
  4 * (pp_2 - pp_1)
 16 * (pp_3 - pp_2)
 64 * (pp_4 - pp_3)

If the results from w1dcn are second-order accurate, then the four curves should (roughly) coincide, with better and better agreement as the resolution increases (i.e. for those curves with a higher density of markers). This is precisely what we see, and this test thus provides a very strong indication that the code is second-order convergent. (In order to establish that convergence is to a solution of the original PDE, we could use the technique of independent residual evaluation; see section 1.9 of the notes on Basic Finite Difference Techniques for an explanation).

IMPORTANT FINAL NOTE!

Due to the fact that RNPL evolutions always start from an initial state file (such as in0.sdf), one must be very careful that the state file that is read by the application actually corresponds to the desired run parameters. It is a common error to run an RNPL program with the wrong initial state file, and, unfortunately, if can sometimes be difficult to detect that this has happened.

You should get in the habit of always removing the initial state file before running an RNPL code, and you may wish to consider defining a shell alias or script to combine the operations of removing the state file and running the RNPL program.