## Interactive session for prepared example 2

```
*******************************************************
****  Welcome to the PDE2D 9.6 Interactive Driver  ****
*******************************************************

PDE2D can solve systems of (linear or nonlinear) steady-state,
time-dependent and eigenvalue partial differential equations
in 1D intervals, general 2D regions, and in a wide range of simple
3D regions.  Ordinary differential equation systems can also be
solved.

You will now be asked a series of interactive questions about your
problem.  The answers you give will be used to construct a PDE2D
FORTRAN program, which can then be compiled and linked with the PDE2D
runtime routines to produce an executable program.  The FORTRAN
driver program created will be well-documented and highly readable
(most of the interactive messages are repeated in the comments),
so that minor modifications or corrections can be made directly to
the FORTRAN program, without the need to work through a new
interactive session.

You can alternatively create your PDE2D FORTRAN program using the
PDE2D graphical user interface (GUI) ("pde2d_gui [progname]").  It
[RETURN]
is extraordinarily easy to set up problems using the PDE2D GUI,
which handles 0D and 1D problems, and 2D and 3D problems in "a wide
range of simple regions".  However, the PDE2D GUI cannot handle
complex regions, so if you have a complex 2D region you must use
this Interactive Driver.

If this is your first time to use PDE2D, you may want to work
through an example problem before trying one of your own.  Do you
want to work through a prepared example?
|---- Enter yes or no
yes
Several prepared examples are available.  Enter:

1 - to see a simple problem: a simply-supported elastic plate
equation, with a unit load concentrated at the midpoint of
a square.
2 - to see a more complex problem: a non-linear, steady-state
PDE, solved in an annulus.  Dirichlet (U = ...) boundary
conditions are imposed on part of the boundary, and Neumann
(dU/dn = ...) conditions are imposed on the other part, in
this example.  The initial triangulation is generated
automatically, and adaptive grid refinement is illustrated.
3 - to see an eigenvalue problem.  The region has a curved
interface across which material properties vary abruptly, in
this example.
4 - to see the first part of a thermal stress problem.  In this
part, the temperature distribution in a V-notched block is
calculated by solving the time-dependent heat conduction
equation, using adaptive time step control.
5 - to see the second part of a thermal stress problem.  In
this part, the stresses induced in the V-notched block by
thermal expansion are calculated, using the temperature
distribution output by example 4.  You must run example 4
[RETURN]
and save the tabular output before you can run example 5.
Examples 4 and 5 illustrate communication between problems.
6 - to see a 1D time-dependent integro-differential equation for
a financial math application.  In this problem there is a
term involving an integral of the solution, which requires
that we use PDE2D's feature for interpolating the solution
at the last saved time step, for use in the integral term.
7 - to see a waveguide problem (an eigenvalue problem in which
the eigenvalue appears nonlinearly).  This example shows how
to handle boundary conditions of different types on the same
arc, and how to produce a plot of a computed integral vs time.
8 - to see the Navier-Stokes equations (penalty formulation)
solved for a fluid flowing around a bend.
9 - to see a 3D elasticity problem, solved in a torus.  This
example illustrates the use of user-defined coordinate
transformations to handle more general 3D regions.
10 - to see a time-dependent wave equation (reduced to a system
of two PDEs), solved in a 3D box.
11 - to see a 3D eigenvalue problem (the Schrodinger equation
in a hydrogen atom).  This example illustrates the use of
spherical coordinates and periodic boundary conditions.
12 - to see a 3D eigenvalue problem, solved in a composite region
[RETURN]
consisting of two cylinders of different material properties.
13 - to see the axisymmetric Navier-Stokes equations solved in a
non-rectangular channel, using the collocation FEM.
14 - to see a 1D saturated/unsaturated water flow problem.
15 - to see a 1D version of the Schrodinger eigenvalue equation
of example 11.
0 - (no example)

If you select one of the examples, the correct answer for each
interactive question will be supplied after the question.
|---- Enter an integer value in the range 0 to 15
2
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ The example problem is                                               \$
\$                                                                      \$
\$      Uxx + Uyy = U**3                                                \$
\$                                                                      \$
\$ in an annulus with inner radius R0=0.05, outer radius R1=1, with     \$
\$ dU/dn = -1/R1**2 on the outer boundary, U=1/R0 on the inner          \$
\$ boundary.                                                            \$
\$                                                                      \$
\$ where Uxx=second derivative of U with respect to X, etc.             \$
\$     dU/dn=normal derivative of U.                                    \$
\$ The exact solution is U = 1/sqrt(X**2+Y**2).                         \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
[RETURN]
In what follows, when you are told to enter a 'FORTRAN expression',
this means any valid FORTRAN expression of 65 characters or less.
In this expression, you may include references to FORTRAN function
subprograms.  You may define these functions line by line at the end
of the interactive session, when prompted, or add them later using an
editor.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you enter a "#" in the first column of any input line, this       +
+ instructs the interactive driver to read this and subsequent input   +
+ lines from the file "pde2d.in".  A "#" in the first column of an     +
+ input line in the file "pde2d.in" (or an end-of-file) instructs the  +
+ driver to switch back to interactive input.                          +
+                                                                      +
+ All lines input during an interactive session are echo printed to    +
+ a file "echo.out".  You may want to modify this file and rename it   +
+ "pde2d.in", and read some or all of your input from this file during +
+ your next interactive session.                                       +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++

PDE2D can solve problems with 0,1,2 or 3 space variables.  Enter the
dimension of your problem:
[RETURN]

0 - to solve a time-dependent ordinary differential equation system,
or an algebraic or algebraic eigenvalue system
1 - to solve problems in 1D intervals
2 - to solve problems in general 2D regions
3 - to solve problems in a wide range of simple 3D regions

\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: 2                                                             \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter an integer value in the range 0 to 3
2
Which finite element method do you want to use:

1. Galerkin method
2. Collocation method

The collocation method can handle a wide range of simple 2D regions;
the Galerkin method can handle completely general 2D regions.

Enter 1 or 2 to select a finite element method.
+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you have a composite region, with discontinuous material          +
+ parameters, you should use the Galerkin method.  If your partial     +
+ differential equations and boundary conditions are difficult to put  +
+ into the "divergence" form required by the Galerkin method, or if    +
+ you have periodic boundary conditions, use the collocation method.   +
+ The collocation method produces an approximate solution with         +
+ continuous first derivatives; the Galerkin solution is continuous    +
+ but its first derivatives are not.                                   +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: 1                                                             \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
[RETURN]
|---- Enter an integer value in the range 1 to 2
1
Is double precision mode to be used?  Double precision is recommended.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If double precision mode is used, variables and functions assigned   +
+ names beginning with a letter in the range A-H or O-Z will be DOUBLE +
+ PRECISION, and you should use double precision constants and FORTRAN +
+ expressions throughout; otherwise such variables and functions will  +
+ be of type REAL.  In either case, variables and functions assigned   +
+ names beginning with I,J,K,L,M or N will be of INTEGER type.         +
+                                                                      +
+ It is possible to convert a single precision PDE2D program to double +
+ precision after it has been created, using an editor.  Just change   +
+ all occurrences of "real" to "double precision"                      +
+                    " tdp" to "dtdp"  (note leading blank)            +
+ Any user-written code or routines must be converted "by hand", of    +
+ course.  To convert from double to single, reverse the changes.      +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
If you don't want to read the FINE PRINT, default NPROB.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you want to solve several similar problems in the same run, set   +
+ NPROB equal to the number of problems you want to solve.  Then NPROB +
+ loops through the main program will be done, with IPROB=1,...,NPROB, +
+ and you can make the problem parameters vary with IPROB.  NPROB      +
+ defaults to 1.                                                       +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ We want to run the problem twice, to do adaptive triangulation       \$
\$ refinement, so                                                       \$
\$ enter: NPROB = 2                                                     \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NPROB =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
2
PDE2D solves the time-dependent system (note: U,A,B,F,FB,GB,U0 may be
vectors, C,RHO may be matrices):

C(X,Y,T,U,Ux,Uy)*d(U)/dT = d/dX* A(X,Y,T,U,Ux,Uy)
+ d/dY* B(X,Y,T,U,Ux,Uy)
-       F(X,Y,T,U,Ux,Uy)

or the steady-state system:

d/dX* A(X,Y,U,Ux,Uy)
+ d/dY* B(X,Y,U,Ux,Uy)
= F(X,Y,U,Ux,Uy)

or the linear and homogeneous eigenvalue system:

d/dX* A(X,Y,U,Ux,Uy)
+ d/dY* B(X,Y,U,Ux,Uy)
= F(X,Y,U,Ux,Uy) + lambda*RHO(X,Y)*U

in an arbitrary two-dimensional region, R, with 'fixed' boundary
conditions on part of the boundary:

[RETURN]
U = FB(X,Y,[T])

and 'free' boundary conditions on the other part:

A*nx + B*ny = GB(X,Y,[T],U,Ux,Uy)

For time-dependent problems there are also initial conditions:

U = U0(X,Y)   at T=T0

Here Ux,Uy represent the (vector) functions dU/dX,dU/dY, and (nx,ny)
represents the unit outward normal to the boundary.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If your PDEs involve the solution at points other than (X,Y), the    +
+ function                                                             +
+             (D)OLDSOL2(IDER,IEQ,XX,YY,KDEG)                          +
+ will interpolate (using interpolation of degree KDEG=1,2 or 3) to    +
+ (XX,YY) the function saved in UOUT(*,*,IDER,IEQ,ISET) on the last    +
+ time step or iteration (ISET) for which it has been saved.  Thus,    +
+ for example, if IDER=1, this will return the latest value of         +
+ component IEQ of the solution at (XX,YY), assuming this has not been +
[RETURN]
+ modified using UPRINT...  If your equations involve integrals of the +
+ solution, for example, you can use (D)OLDSOL2 to approximate these   +
+ using the solution from the last time step or iteration.             +
+                                                                      +
+ CAUTION: For a steady-state or eigenvalue problem, you must reset    +
+ NOUT=1 if you want to save the solution each iteration.              +
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ A system of NEQN complex partial differential equations must be      +
+ written as a system of 2*NEQN real equations, by separating the      +
+ equations into their real and imaginary parts.  However, note that   +
+ the complex arithmetic abilities of FORTRAN can be used to simplify  +
+ this separation.  For example, the complex PDE:                      +
+       I*(Uxx+Uyy) = 1/(1+U**10),  where U = UR + UI*I                +
+ would be difficult to split up analytically, but using FORTRAN       +
+ expressions it is easy:                                              +
+   A1 = -UIx,  B1 = -UIy,  F1 =  REAL(1.0/(1.0+CMPLX(UR,UI)**10))     +
+   A2 =  URx,  B2 =  URy,  F2 = AIMAG(1.0/(1.0+CMPLX(UR,UI)**10))     +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
You may now define global parameters, which may be referenced in any
of the "FORTRAN expressions" you input throughout the rest of this
interactive session.  You will be prompted alternately for parameter
names and their values; enter a blank name when you are finished.
[RETURN]

Parameter names are valid FORTRAN variable names, starting in
column 1.  Thus each name consists of 1 to 6 alphanumeric characters,
the first of which must be a letter.  If the first letter is in the
range I-N, the parameter must be an integer.

Parameter values are either FORTRAN constants or FORTRAN expressions
involving only constants and global parameters defined on earlier
lines.  They may also be functions of the problem number IPROB, if
you are solving several similar problems in one run (NPROB > 1).  Note
that you are defining global CONSTANTS, not functions; e.g., parameter
values may not reference any of the independent or dependent variables

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you define other parameters here later, using an editor, you must +
+ add them to COMMON block /PARM8Z/ everywhere this block appears, if  +
+ they are to be "global" parameters.                                  +
+                                                                      +
+ The variable PI is already included as a global parameter, with an   +
+ accurate value 3.14159...                                            +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
[RETURN]
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: R0                                                            \$
\$          0.05                                                        \$
\$        R1                                                            \$
\$          1                                                           \$
\$      [blank line]                                                    \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
Parameter name = (type blank line to terminate)
R0
R0     =
|----Enter constant or FORTRAN expression-----------------------|
0.05
Parameter name = (type blank line to terminate)
R1
R1     =
|----Enter constant or FORTRAN expression-----------------------|
1
Parameter name = (type blank line to terminate)

If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Do you want to be given a chance to write a FORTRAN block before the +
+ definitions of many functions?  If you answer 'no', you will still   +
+ be given a chance to write code before the definition of the PDE     +
+ coefficients, but not other functions.  Of course, you can always    +
+ add code later directly to the resulting program, using an editor.   +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
Now, the first step in creating a 2D PDE2D program is to construct
an initial triangulation of the region R over which the partial
differential equations are to be solved.  This initial triangulation
can later be refined and graded to your specifications.

The boundary of the region R should be divided into distinct
(curved or straight) arcs, each of which is smooth with smooth
boundary conditions.  Thus at every corner or point where the boundary
conditions have a discontinuity or change type, a new boundary arc
should begin.  Each arc is assigned a unique integer arc number, which
is arbitrary except that it must be negative if 'fixed' boundary
conditions are specified on that arc, and positive (but < 1000) if
'free' boundary conditions are specified.

'Fixed' means that all unknowns are specified on that boundary
(Dirichlet type conditions) and 'free' refers to more general boundary
conditions, as defined earlier.  These more general conditions include
specification of the boundary flux (or normal derivative in the case
of a Laplacian operator).

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If some, but not all, of the unknowns are specified on a boundary    +
[RETURN]
+ arc, that arc should be considered 'free', since even 'fixed'        +
+ type conditions of the form:                                         +
+                   Ui = FBi(X,Y,[T])                                  +
+ can be expressed as 'free' boundary conditions in the form:          +
+               Ai*nx + Bi*ny = zero(Ui-FBi(X,Y,[T]))                  +
+ where zero(f) = Big_Number*f is a PDE2D-supplied function.           +
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ If Ai=Bi=0 (eg, if the PDE is first order), then setting GBi=0 is    +
+ equivalent to setting "no" boundary condition (which is sometimes    +
+ appropriate for first order PDEs), because the boundary condition    +
+ Ai*nx + Bi*ny = GBi reduces to 0=0.                                  +
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Note: while 'fixed' boundary conditions are enforced exactly (at the +
+ nodes), 'free' boundary conditions are not; the greater the overall  +
+ solution accuracy, the more closely they are satisfied.              +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++

There are 3 options for generating the initial triangulation:

1. If the region R is a rectangle with sides parallel to the X and Y
axes, then the initial triangulation can be generated automatically.
2. If the region R can be conveniently described by parametric
[RETURN]
equations in the form:
X = X(P,Q)
Y = Y(P,Q)
where P and Q have constant limits (for example, if X=P*COS(Q),
Y=P*SIN(Q), P and Q will represent polar coordinates), then the
initial triangulation can also be generated automatically.
3. For more general regions, you will need to create an initial
triangulation by hand.

Enter INTRI = 1,2 or 3 to choose an initial triangulation option.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ We can parameterize using polar coordinates, with an arc number of   \$
\$ -1 assigned to the inner boundary, +1 to the outer boundary.         \$
\$                                                                      \$
\$ enter: INTRI = 2                                                     \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
INTRI =
|---- Enter an integer value in the range 1 to 3
2
Define the Cartesian coordinates X,Y as functions of the parameters
P and Q, so that your region is generated as P and Q vary between
constant limits PGRID(1) < P < PGRID(NPGRID)
QGRID(1) < Q < QGRID(NQGRID)
where PGRID and QGRID will be defined later.  X=P, Y=Q by default.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ X and Y must be continuous functions of P and Q, but their           +
+ derivatives need not be continuous.                                  +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ The region can be described using polar coordinates, with radius P   \$
\$ and angle Q, with R0 < P < R1 and 0 < Q < 2*PI.                      \$
\$ enter: X = P*COS(Q)                                                  \$
\$        Y = P*SIN(Q)                                                  \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
X =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
P*COS(Q)
Y =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
P*SIN(Q)
For a parameterized region, the initial triangulation is defined by a
set of curved grid lines corresponding to
P = PGRID(1),PGRID(2),...,PGRID(NPGRID)
Q = QGRID(1),QGRID(2),...,QGRID(NQGRID)
Each grid quadrilateral is divided into three or four triangles.

You will first be prompted for NPGRID, the number of P-grid points,
then for PGRID(1),...,PGRID(NPGRID).  Any points defaulted will be
uniformly spaced between the points you define; the first and last
represent the lower and upper limits of P, and cannot be defaulted.
Then you will be prompted similarly for the number and values of
the Q-grid points.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: NPGRID = 4                                                    \$
\$        PGRID(1) = R0                                                 \$
\$        PGRID(NPGRID) = R1                                            \$
\$ and default PGRID(2),PGRID(3)                                        \$
\$        NQGRID = 9                                                    \$
\$        QGRID(1) = 0.0                                                \$
\$        QGRID(NQGRID) = 2*PI                                          \$
\$ and default QGRID(2),...,QGRID(8)                                    \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
[RETURN]
NPGRID =
|---- Enter an integer value in the range 2 to +INFINITY
4
PGRID(1) =
|----Enter constant or FORTRAN expression-----------------------|
R0
PGRID(2) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

PGRID(3) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

PGRID(NPGRID) =
|----Enter constant or FORTRAN expression-----------------------|
R1
NQGRID =
|---- Enter an integer value in the range 2 to +INFINITY
9
QGRID(1) =
|----Enter constant or FORTRAN expression-----------------------|
0.0
QGRID(2) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

QGRID(3) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

QGRID(4) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

QGRID(5) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

QGRID(6) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

QGRID(7) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

QGRID(8) =              (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

QGRID(NQGRID) =
|----Enter constant or FORTRAN expression-----------------------|
2*PI
Enter the arc numbers IPARC(1),IPARC(2) corresponding to the arcs
P = PGRID(1) and P = PGRID(NPGRID), respectively, and the arc numbers
IQARC(1),IQARC(2) corresponding to the arcs Q = QGRID(1) and
Q = QGRID(NQGRID).  Recall that negative arc numbers correspond to
'fixed' boundary conditions, and positive arc numbers ( < 1000)
correspond to 'free' boundary conditions.

For a parameterized region, an "arc" may be a single point (for
example, R=0 if polar coordinates are used), in which case the arc
number must be set to 0, or it may be interior to the region (for
example, Theta=0 and Theta=2*pi if polar coordinates are used), in
which case the arc number must be > 999.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: IPARC(1) = -1     IPARC(2) = 1                                \$
\$        IQARC(1) = 1001   IQARC(2) = 1001                             \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
IPARC(1) =
|----Enter constant or FORTRAN expression-----------------------|
-1
IPARC(2) =
|----Enter constant or FORTRAN expression-----------------------|
1
IQARC(1) =
|----Enter constant or FORTRAN expression-----------------------|
1001
IQARC(2) =
|----Enter constant or FORTRAN expression-----------------------|
1001
How many triangles (NTF) are desired for the final triangulation?
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: NTF = 1000                                                    \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NTF =
|---- Enter an integer value in the range 96 to +INFINITY
1000
Enter a FORTRAN expression for TRIDEN(X,Y), which controls the
grading of the triangulation.  TRIDEN should be largest where the
triangulation is to be most dense.  The default is TRIDEN(X,Y)=1.0
(a uniform triangulation).

TRIDEN may also be a function of the initial triangle number KTRI.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ press [RETURN] to default TRIDEN                                     \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
TRIDEN =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Do you want the triangulation to be graded adaptively?               +
+                                                                      +
+ If you answer "yes", make sure there is no "pde2d.adp" file in the   +
+ working directory the first time you run the program, then run the   +
+ program two or more times, possibly increasing NTF each time.  On    +
+ the first run, the triangulation will be graded as guided by         +
+ TRIDEN(X,Y), but on each subsequent run, information output by the   +
+ previous run to "pde2d.adp" will be used to guide the grading of the +
+ new triangulation.                                                   +
+
+ If ADAPT=.TRUE., after each run a file "pde2d.adp" is written        +
+ which tabulates the values of the magnitude of the gradient of the   +
+ solution (at the last time step or iteration) at an output NXP8Z by  +
+ NYP8Z grid of points (NXP8Z and NYP8Z are set to 101 in a PARAMETER  +
+ statement in the main program, so they can be changed if desired).   +
+ If NEQN > 1, a normalized average of the gradients of the NEQN       +
+ solution components is used.                                         +
+                                                                      +
+ You can do all the "runs" in one program, by setting NPROB > 1.      +
[RETURN]
+ Each pass through the DO loop, PDE2D will read the gradient values   +
+ output the previous pass.  If RESTRT=.TRUE., GRIDID=.FALSE., and     +
+ T0,TF are incremented each pass through the DO loop, it is possible  +
+ in this way to solve a time-dependent problem with an adaptive,      +
+ moving, grid.                                                        +
+                                                                      +
+ Increase the variable EXAG from its default value of 1.5 if you want +
+ to exaggerate the grading of an adaptive triangulation (make it less +
+ uniform).  EXAG should normally not be larger than about 2.0.        +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++

If you don't want to read the FINE PRINT, enter 'no'.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
If you don't want to read the FINE PRINT, default SHAPE.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Enter a FORTRAN expression for SHAPE(X,Y), which controls the        +
+ approximate shape of the triangles.  The triangulation refinement    +
+ will proceed with the goal of generating triangles with an average   +
+ height to width ratio of approximately SHAPE(X,Y) near the point     +
+ (X,Y).  SHAPE must be positive.  The default is SHAPE(X,Y)=1.0.      +
+                                                                      +
+ SHAPE may also be a function of the initial triangle number KTRI.    +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ press [RETURN] to default SHAPE                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
SHAPE =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

If you don't want to read the FINE PRINT, enter ISOLVE = 4.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ The following linear system solvers are available:                   +
+                                                                      +
+ 1. Band method                                                       +
+               The band solver uses a reverse Cuthill-McKee ordering. +
+ 2. Frontal method                                                    +
+               This is an out-of-core version of the band solver.     +
+ 3. Jacobi bi-conjugate gradient method                               +
+               This is a preconditioned bi-conjugate gradient, or     +
+               Lanczos, iterative method.  (This solver is MPI-       +
+               enhanced, if MPI is available.)  If you want to        +
+               override the default convergence tolerance, set a      +
+               new relative tolerance CGTL8Z in the main program.     +
+ 4. Sparse direct method                                              +
+               This is based on Harwell Library routines MA27/MA37,   +
+               developed by AEA Industrial Technology at Harwell      +
+               Laboratory, Oxfordshire, OX11 0RA, United Kingdom      +
+               (used by permission).                                  +
+ 5. Local solver                                                      +
+               Choose this option ONLY if alternative linear system   +
[RETURN]
+               solvers have been installed locally.  See subroutines  +
+               (D)TD3M, (D)TD3N in file (d)subs.f for instructions    +
+               on how to add local solvers.                           +
+ 6. MPI-based parallel band solver                                    +
+               This is a parallel solver which runs efficiently on    +
+               multiple processor machines, under MPI.  It is a       +
+               band solver, with the matrix distributed over the      +
+               available processors.  Choose this option ONLY if the  +
+               solver has been activated locally.  See subroutine     +
+               (D)TD3O in file (d)subs.f for instructions on how to   +
+               activate this solver and the MPI-enhancements to the   +
+               conjugate gradient solver.                             +
+                                                                      +
+ Enter ISOLVE = 1,2,3,4,5 or 6 to select a linear system solver.      +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++

If you don't want to read the FINE PRINT, enter ISOLVE = 4.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: ISOLVE = 4                                                    \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
ISOLVE =
|---- Enter an integer value in the range 1 to 6
4
Enter the element degree (1,2,3 or 4) desired.  A suggested value is
IDEG = 3.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ A negative value for IDEG can be entered, and elements of degree     +
+ ABS(IDEG) will be used, with a lower order numerical integration     +
+ scheme.  This results in a slight increase in speed, but negative    +
+ values of IDEG are normally not recommended.                         +
+                                                                      +
+ The spatial discretization error is O(h**2), O(h**3), O(h**4) or     +
+ O(h**5) when IDEG = 1,2,3 or 4, respectively, is used, where h is    +
+ the maximum triangle diameter, even if the region is curved.         +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: IDEG = 3                                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
IDEG =
|---- Enter an integer value in the range -4 to 4
3
What type of PDE problem do you want to solve?

1. a steady-state (time-independent) problem
2. a time-dependent problem
3. a linear, homogeneous eigenvalue problem

Enter 1,2 or 3 to select a problem type.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: 1                                                             \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter an integer value in the range 1 to 3
1
Is this a linear problem? ("linear" means all differential equations
and all boundary conditions are linear).  If you aren't sure, it is
safer to answer "no".
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
Give an upper limit on the number of Newton's method iterations
(NSTEPS) to be allowed for this nonlinear problem.  NSTEPS defaults
to 15.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ The iteration will stop if convergence occurs before the upper       +
+ limit has been reached.  In the FORTRAN program created by the       +
+ preprocessor, NCON8Z will be .TRUE. if Newton's method has converged.+
+                                                                      +
+ For highly non-linear problems you may want to construct a one-      +
+ parameter family of problems using the variable T, such that for     +
+ T=1 the problems is easy (e.g. linear) and for T > N (N is less      +
+ NSTEPS), the problem reduces to the original highly nonlinear        +
+ problem.  For example, the nonlinear term(s) may be multiplied by    +
+ MIN(1.0,(T-1.0)/N).                                                  +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ press [RETURN] to default NSTEPS                                     \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NSTEPS =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

How many differential equations (NEQN) are there in your problem?
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: NEQN = 1                                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NEQN =
|---- Enter an integer value in the range 1 to 99
1
You may now choose names for the component(s) of the (possibly vector)
solution U.  Each must be an alphanumeric string of one to three
characters, beginning with a letter in the range A-H or O-Z.  The
variable names X,Y,T,S,A and B must not be used.  The name should
start in column 1.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: U1 = U                                                        \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
U1 =
U
PDE2D solves the equation:

d/dX* A(X,Y,U,Ux,Uy)
+ d/dY* B(X,Y,U,Ux,Uy)
= F(X,Y,U,Ux,Uy)

with 'fixed' boundary condition:

U = FB(X,Y)

or 'free' boundary condition:

A*nx + B*ny = GB(X,Y,U,Ux,Uy)

where U(X,Y) is the unknown and F,A,B,FB,GB are user-supplied
functions.

note:
(nx,ny) = unit outward normal to the boundary
Ux = d(U)/dX
Uy = d(U)/dY

[RETURN]
Is this problem symmetric?  If you don't want to read the FINE PRINT,
it is safe to enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ This problem is called symmetric if the matrix                       +
+                                                                      +
+                F.U     F.Ux     F.Uy                                 +
+                A.U     A.Ux     A.Uy                                 +
+                B.U     B.Ux     B.Uy                                 +
+                                                                      +
+ is always symmetric, where F.U means d(F)/d(U), and similarly        +
+ for the other terms.  In addition, GB must not depend on Ux,Uy.      +
+                                                                      +
+ The memory and execution time are halved if the problem is known to  +
+ be symmetric.                                                        +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ The Jacobian matrix is diagonal so,                                  \$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
If you don't want to read the FINE PRINT, enter 'yes' (strongly
recommended).

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ The partial derivatives of some of the PDE and boundary condition    +
+ coefficients are required by PDE2D.  These may be calculated         +
+ automatically using a finite difference approximation, or supplied   +
+ by the user.  Do you want them to be calculated automatically?       +
+                                                                      +
+ If you answer 'yes', you will not be asked to supply the derivatives,+
+ but there is a small risk that the inaccuracies introduced by the    +
+ finite difference approximation may cause the Newton iteration       +
+ to converge more slowly or to diverge, especially if low precision   +
+ is used.  This risk is very low, however, and since answering 'no'   +
+ means you may have to compute many partial derivatives, it is        +
+ recommended you answer 'yes' unless you have some reason to believe  +
+ there is a problem with the finite difference approximations.        +
+                                                                      +
+ If you supply analytic partial derivatives, PDE2D will do some spot  +
+ checking and can usually issue a warning if any are supplied         +
+ incorrectly.                                                         +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
[RETURN]
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
You may calculate one or more integrals (over the entire region) of
some functions of the solution and its derivatives.  How many integrals
(NINT), if any, do you want to calculate?

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ In the FORTRAN program created by the preprocessor, the computed     +
+ values of the integrals will be returned in the vector SINT8Z.  If   +
+ several iterations or time steps are done, only the last computed    +
+ values are saved in SINT8Z (all values are printed).                 +
+                                                                      +
+ A limiting value, SLIM8Z(I), for the I-th integral can be set        +
+ below in the main program.  The computations will then stop          +
+ gracefully whenever SINT8Z(I) > SLIM8Z(I), for any I=1...NINT.       +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: NINT = 1                                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NINT =
|---- Enter an integer value in the range 0 to 20
1
Enter FORTRAN expressions for the functions whose integrals are to be
calculated and printed.  They may be functions of

X,Y,U,Ux,Uy and (if applicable) T

The integrals may also contain references to the initial triangle
number KTRI.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you only want to integrate a function over part of the region,    +
+ define that function to be zero in the rest of the region.           +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ We want to calculate the L1 norm of the error, so                    \$
\$ enter: INTEGRAL = ABS(U-TRUE(X,Y))                                   \$
\$   where TRUE(X,Y) is the true solution, to be defined later.         \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
INTEGRAL =
|----Enter constant or FORTRAN expression-----------------------|
ABS(U-TRUE(X,Y))
You may calculate one or more boundary integrals (over the entire
boundary) of some functions of the solution and its derivatives.  How
many boundary integrals (NBINT), if any, do you want to calculate?

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ In the FORTRAN program created by the preprocessor, the computed     +
+ values of the integrals will be returned in the vector BINT8Z.  If   +
+ several iterations or time steps are done, only the last computed    +
+ values are saved in BINT8Z (all values are printed).                 +
+                                                                      +
+ A limiting value, BLIM8Z(I), for the I-th boundary integral can be   +
+ set below in the main program.  The computations will then stop      +
+ gracefully whenever BINT8Z(I) > BLIM8Z(I), for any I=1...NBINT.      +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: NBINT = 0                                                     \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NBINT =
|---- Enter an integer value in the range 0 to 20
0
If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Do you have point source terms in your PDEs, involving "Dirac Delta" +
+ functions?  A Dirac Delta function DEL(x-xd0,y-yd0) is a function    +
+ whose integral is 1, but which is zero everywhere except at a        +
+ singular point (xd0,yd0).                                            +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
Now enter FORTRAN expressions to define the PDE coefficients, which
may be functions of

X,Y,U,Ux,Uy

They may also be functions of the initial triangle number KTRI
and, in some cases, of the parameter T.

Recall that the PDE has the form

d/dX*A + d/dY*B = F

\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ The partial differential equation may be written                     \$
\$                                                                      \$
\$  d/dX*(Ux) + d/dY*(Uy) = U**3                                        \$
\$                                                                      \$
\$ When asked if you want to write a FORTRAN block,                     \$
\$ enter: no                                                            \$
\$   then enter the following, when prompted:                           \$
\$    F = U**3         A = Ux         B = Uy                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
[RETURN]

Do you want to write a FORTRAN block to define some parameters to be
used in the definition of these coefficients?
|---- Enter yes or no
no
F =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
U**3
A =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
Ux
B =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
Uy
Now the initial values for Newton's method must be defined using
FORTRAN expressions.  They may be functions of X and Y, and of the
initial triangle number KTRI.

It is important to provide initial values which are at least of the
correct order of magnitude.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: U0 = 1.0                                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
U0 =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
1.0
If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Do you want to read the initial conditions from the restart file,    +
+ if it exists (and use the conditions supplied above if it does not   +
+ exist)?                                                              +
+                                                                      +
+ If so, PDE2D will dump the final solution at the end of each run     +
+ into a restart file "pde2d.res".  Thus the usual procedure for       +
+ using this dump/restart option is to make sure there is no restart   +
+ file in your directory left over from a previous job, then the       +
+ first time you run this job, the initial conditions supplied above   +
+ will be used, but on the second and subsequent runs the restart file +
+ from the previous run will be used to define the initial conditions. +
+                                                                      +
+ You can do all the "runs" in one program, by setting NPROB > 1.      +
+ Each pass through the DO loop, T0,TF,NSTEPS and possibly other       +
+ parameters may be varied, by making them functions of IPROB.         +
+                                                                      +
+ If the 2D or 3D collocation method is used, the coordinate           +
+ transformation should not change between dump and restart.           +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
[RETURN]
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ We want to restart the second solution (IPROB=2) from the first      \$
\$ solution, so                                                         \$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
If the program dumping the solution and the program restarting from
the dump have identical finite element grids, the solution at the
nodes can be dumped to the restart file "pde2d.res", and there will
be no loss of accuracy due to interpolation when this solution is read
back.  Will the dumping and restarting programs have identical grids?

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you answer 'no', then you can change the grid between runs, but   +
+ the solution on the new grid will be interpolated from the dumped    +
+ solution, so there may be some loss of accuracy due to interpolation.+
+ The solution will be dumped to the restart file on a uniform grid    +
+ of size                                                              +
+    for 1D problems:   NXP8Z=1001                                     +
+        2D problems:   NXP8Z=101 by NYP8Z=101                         +
+ The values of NXP8Z,NYP8Z are set in a PARAMETER statement in the    +
+ main program, so they can be changed, if desired.  The interpolation +
+ degree (KDEG8Z) is also set in this PARAMETER statement, to 1        +
+ (linear interpolation), but KDEG8Z can be reset to 2 or 3 if desired.+
+                                                                      +
+ If you want to modify the data in the restart file after it has      +
+ been dumped (each trip through the loop IPROB=1,NPROB in the main    +
+ program), answer 'no' even if your grid does not change, then        +
[RETURN]
+ uncomment the two calls to (D)TDPR(1,2) in the main program (and     +
+ also the dimension statement for xres8z,...,ures8z).  You can modify +
+ the restart file by modifying the solution as stored in the array    +
+ URES8Z, after it is read from 'pde2d.res' by the first call, and     +
+ before it is written back by the second call.  The solution is       +
+ stored as follows:                                                   +
+          (D)TDPR1:  URES8Z(IEQN,I) = unknown #IEQN at                +
+                                X = XRES8Z(I)                         +
+          (D)TDPR2:  URES8Z(IEQN,I,J) = unknown #IEQN at              +
+                                X = XRES8Z(I), Y = YRES8Z(J))         +
+ where I=1,...,NXP8Z, J=1,...,NYP8Z.  You may NOT modify the arrays   +
+ XRES8Z,YRES8Z.                                                       +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++

Will the dumping and restarting programs have identical grids?
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ Since adaptive triangulation is requested, the triangulation may     \$
\$ change between IPROB=1 and IPROB=2, so                               \$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
Now the 'fixed' boundary conditions are described.

Are there any boundary arcs with negative arc numbers, where nonzero
'fixed' boundary conditions must be specified?  (The boundary
conditions on arcs with negative numbers default to FB=0.)
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
Enter the arc number (IARC) of a boundary arc with negative arc number.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: IARC = -1                                                     \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
IARC =
|---- Enter an integer value in the range -INFINITY to -1
-1
Enter a FORTRAN expression to define FB on this arc.  It may be a
function of X,Y and (if applicable) T.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ This function may also reference the initial triangle number KTRI,   +
+ and the arc parameter S (S = P or Q when INTRI=2).                   +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++

Recall that fixed boundary conditions have the form

U = FB

\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: FB = 1.0/R0                                                   \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
FB =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
1.0/R0
Are there any more negative arcs, with nonzero boundary conditions?
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
Now the 'free' boundary conditions are described.

Are there any boundary arcs with positive (but < 1000) arc numbers,
where nonzero 'free' boundary conditions must be specified?  (The
boundary conditions on arcs with positive numbers default to GB=0.)
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
Enter the arc number (IARC) of a boundary arc with positive arc number.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: IARC = 1                                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
IARC =
|---- Enter an integer value in the range 1 to 999
1
Enter FORTRAN expressions to define the following free boundary
condition functions on this arc.  They may be functions of

X,Y,U,Ux,Uy and (if applicable) T

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ These functions may also reference the components (NORMx,NORMy) of   +
+ the unit outward normal vector, the initial triangle number KTRI,    +
+ and the arc parameter S (S = P or Q when INTRI=2).                   +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++

Recall that free boundary conditions have the form

A*nx+B*ny = GB

\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ Recall that GB = A*nx+B*ny = Ux*nx+Uy*ny = dU/dn = -1/R1**2          \$
\$ Thus enter, when prompted:                                           \$
\$    GB = -1.0/R1**2                                                   \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
GB =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
-1.0/R1**2
Are there any more positive arcs, with nonzero boundary conditions?
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
If you don't want to read the FINE PRINT, default all of the following
variables.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Normally, PDE2D saves the values of U,A,B at the output points.      +
+ If different variables are to be saved (for later printing or        +
+ plotting) the following functions can be used to re-define the       +
+ output variables:                                                    +
+    define UPRINT(1) to replace U                                     +
+           APRINT(1)            A                                     +
+           BPRINT(1)            B                                     +
+ Each function may be a function of                                   +
+                                                                      +
+           X,Y,U,Ux,Uy,A,B and (if applicable) T                      +
+                                                                      +
+ Each may also be a function of the initial triangle number KTRI and  +
+ the integral estimates SINT(1),...,BINT(1),...                       +
+                                                                      +
+ The default for each variable is no change, for example, UPRINT(1)   +
+ defaults to U.  Enter FORTRAN expressions for each of the            +
+ following functions (or default).                                    +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
[RETURN]
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ press [RETURN] to default all output modification variables          \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
Replace U for postprocessing?
UPRINT(1) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

Replace A for postprocessing?
APRINT(1) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

Replace B for postprocessing?
BPRINT(1) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

The solution is normally saved on a NX+1 by NY+1 rectangular grid of
points
(XA + I*(XB-XA)/NX , YA + J*(YB-YA)/NY)
I=0,...,NX, J=0,...,NY.  Enter values for NX and NY.  Suggested values
are NX=NY=25.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you want to save the solution at an arbitrary user-specified      +
+ set of points, set NY=0 and NX+1=number of points.  In this case you +
+ can request tabular output of the solution, but you cannot make any  +
+ solution plots.                                                      +
+                                                                      +
+ If you set NEAR8Z=1 in the main program, the values saved at each    +
+ output point will actually be the solution as evaluated at a nearby  +
+ integration point.  For most problems this obviously will produce    +
+ less accurate output or plots, but for certain (rare) problems, a    +
+ solution component may be much less noisy when plotted only at       +
+ integration points.                                                  +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: NX = 40                                                       \$
\$        NY = 40                                                       \$
[RETURN]
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
NX =
|---- Enter an integer value in the range 1 to +INFINITY
40
NY =
|---- Enter an integer value in the range 0 to +INFINITY
40
The solution is saved on an NX+1 by NY+1 rectangular grid covering the
rectangle (XA,XB) x (YA,YB).  Enter values for XA,XB,YA,YB.  These
variables are usually defaulted.

The default is a rectangle which just covers the entire region.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ press [RETURN] to default XA,XB,YA,YB                                \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
XA =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

XB =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

YA =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

YB =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

To generate tabular or graphical output, choose an output option from
the list below.

0. No further output is desired
1. Table of values at output points
The tabulated output is saved in a file.
2. Surface plot of a scalar variable
3. Contour plot of a scalar variable
4. Vector field plot, with arrows indicating magnitude and direction
5. One-dimensional cross-sectional plots (versus X, Y or T)
or, if applicable:
6. Stress field plot (requires 2 or more PDES)
A plot of the principal stresses. Compression is indicated
by arrows pointing toward each other, and tension by
arrows pointing away from each other.

Enter 0,1,2,3,4,5 or 6 to select an output option.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you set NY=0 and saved the solution at an arbitrary set of        +
+ user-specified output points, you can only request tabular output.   +
+                                                                      +
[RETURN]
+ If you decide later that you want additional types of plots not      +
+ requested during this interactive session, you will have to work     +
+ through a new interactive session, so it is recommended that you     +
+ request all output or plots you think you MIGHT eventually want now, +
+ during this session.                                                 +
+                                                                      +
+ Regardless of the options you select, a dummy subroutine POSTPR      +
+ will be included in the program created by the interactive driver;   +
+ you can add your own postprocessing code to this subroutine.         +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: 3, the first time you see this message and                    \$
\$        0, the second time                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter an integer value in the range 0 to 5
3
Enter a value for IVAR, to select the variable to be plotted or
printed:
IVAR = 1 means U  (possibly as modified by UPRINT,..)
2       A
3       B
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ We want a contour plot of U, so                                      \$
\$ enter: IVAR = 1                                                      \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
IVAR =
|---- Enter an integer value in the range 1 to 3
1
If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Do you want to scale the axes on the plot so that the region is      +
+ undistorted?  Otherwise the axes will be scaled so that the figure   +
+ approximately fills the plot space.                                  +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
Enter lower (UMIN) and upper (UMAX) bounds for the contour values. UMIN
and UMAX are often defaulted.

Labeled contours will be drawn corresponding to the values

UMIN + S*(UMAX-UMIN),    for S=0.05,0.15,...0.95.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ By default, UMIN and UMAX are set to the minimum and maximum values  +
+ of the variable to be plotted.  For a common scaling, you may want   +
+ to set UMIN=ALOW, UMAX=AHIGH.  ALOW and AHIGH are the minimum and    +
+ maximum values over all output points and over all saved time steps  +
+ or iterations.                                                       +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ press [RETURN] to default UMIN and UMAX                              \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
UMIN =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

UMAX =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|

Do you want two additional unlabeled contours to be drawn between each
pair of labeled contours?
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
Enter a title, WITHOUT quotation marks.  A maximum of 40 characters
are allowed.  The default is no title.
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: Solution of second example                                    \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
TITLE =             (Press [RETURN] to default)
|----Enter title or name---------------|
Solution of second example
To generate tabular or graphical output, choose an output option from
the list below.

0. No further output is desired
1. Table of values at output points
The tabulated output is saved in a file.
2. Surface plot of a scalar variable
3. Contour plot of a scalar variable
4. Vector field plot, with arrows indicating magnitude and direction
5. One-dimensional cross-sectional plots (versus X, Y or T)
or, if applicable:
6. Stress field plot (requires 2 or more PDES)
A plot of the principal stresses. Compression is indicated
by arrows pointing toward each other, and tension by
arrows pointing away from each other.

Enter 0,1,2,3,4,5 or 6 to select an output option.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you set NY=0 and saved the solution at an arbitrary set of        +
+ user-specified output points, you can only request tabular output.   +
+                                                                      +
[RETURN]
+ If you decide later that you want additional types of plots not      +
+ requested during this interactive session, you will have to work     +
+ through a new interactive session, so it is recommended that you     +
+ request all output or plots you think you MIGHT eventually want now, +
+ during this session.                                                 +
+                                                                      +
+ Regardless of the options you select, a dummy subroutine POSTPR      +
+ will be included in the program created by the interactive driver;   +
+ you can add your own postprocessing code to this subroutine.         +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: 3, the first time you see this message and                    \$
\$        0, the second time                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter an integer value in the range 0 to 5
0
Do you want to define any FORTRAN function subprograms used in any of
the FORTRAN 'expressions' entered earlier, entering them line by line?

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ If you selected double precision accuracy earlier, be sure to        +
+ declare these functions and their arguments DOUBLE PRECISION.        +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: yes                                                           \$
\$ Then enter one line at a time of the FORTRAN function TRUE:          \$
\$       FUNCTION TRUE(X,Y)                                             \$
\$       IMPLICIT DOUBLE PRECISION (A-H,O-Z)                            \$
\$       TRUE = 1.0/SQRT(X*X+Y*Y)                                       \$
\$       RETURN                                                         \$
\$       END                                                            \$
\$     [blank line]                                                     \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
yes
Remember to begin FORTRAN statements in column 7
|-----7-----Input FORTRAN now (type blank line to terminate)-----------|
FUNCTION TRUE(X,Y)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
TRUE = 1.0/SQRT(X*X+Y*Y)
RETURN
END

If you don't want to read the FINE PRINT, enter 'no'.

+++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
+ Do you want to define any FORTRAN function subprograms used in any   +
+ of the FORTRAN 'expressions' entered earlier, by interpolating the   +
+ tabular output saved in a file created on an earlier PDE2D run?      +
++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
\$ enter: no                                                            \$
\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$ EXAMPLE 2 \$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
|---- Enter yes or no
no
More detailed information about PDE2D can be found in the book
"Solving Partial Differential Equation Applications with PDE2D,"
Granville Sewell, John Wiley & Sons, 2018.

*******************************************
***** Input program has been created  *****
*******************************************

```