NBS

PUBLICATIONS

AlllDE

NBSIR 86-341

NATL INST OF STANDARDS & TECH R.I.C.

A1 11 02247894

Blue, James L/B2DE : a program for solvl QC100 .U56 NO.86-3411 1986 V19 C.1 NBS-P B2DE -- A Program for Solving

Systems of Partial Differential

Equations in Two Dimensions

James L. Blue

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standard’s

Center for Applied Mathematics

Scientific Computing Division

Gaithersburg, MD 20899

June 1 986

100

.056

86-3^11

1986

C. 2

>ARTMENT OF COMMERCE

BUREAU OF STANDARDS

RBS

RESEARCH

INFORMATION CENTER

NBSIR 86-341

B2DE - A PROGRAM FOR SOLVING

SYSTEMS OF PARTIAL DIFFERENTIAL

EQUATIONS IN TWO DIMENSIONS

James L. Blue

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Applied Mathematics

Scientific Computing Division

Gaithersburg, MD 20899

July 1 986

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige, Secretary

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director
mM0flAJ?38 WJTAWBOWJ jJv'M^S-

'^O Bilvl

w

" ¥ I (•<h8 J.i

#'"

?'‘'

i'S^V %.». B''>.

.p‘ "ii'

hi»

' ‘L 1'

ti’.'

V? t1

irl VH’>r>«fP lib t i!i|

tKi;-.

Ira.;

1^.

i---

r f

«>-.J

«nJ(M*«|yl =-0 7M JMT«A*1»a * U i

toMinaiMtv <14? v>*)»u* <^<:<n*4'

I
B2DE - A Program for Solving Systems of Nonlinear Elliptic

Partial Differential Equations in Two Dimensions

Manual and Users’ Guide

James L. Blue

Scientific Computing Division

National Bureau of Standards

Gaithersburg, MD 20899

ABSTRACT

B2DE is a program for solving systems of nonlinear elliptic partial differential

equations (PDEs) in two dimensions. The program is a collection of modules with

an interactive driver. Many types of interactive graphics plots are included.

Users may modify the driver, and may be able to construct a “black box”

program for a restricted class of PDEs.

B2DE is available from the author at the above address.

1. Introduction

B2DE is a program for solving systems of nonlinear elliptic partial differential equations

(PDEs) in two dimensions. B2DE has been used for solving problems as simple as Laplace’s

equation [Kell86| and as complicated as the coupled equations of electrostatics and chaxgedparticle transport in a transistor [Blue83, Wils85].

Since there is no algorithm for solving all solvable nonlinear systems of PDEs, B2DE cannot

work as a “black box” program, as does, for example, a Gaussian elimination program.

Instead, B2DE is a collection of modules; the user must decide what to do at each point

of the solution process. An interactive driver is provided for the user. Expert users may

wish to modify the driver, and may construct a “black box” program for a restricted class

of PDEs. Less expert users may wish to modify some routines designed for easy user

modification.

Several kinds of graphic displays are provided to aid the user in imderstanding the solution

and the solution process.

This manual includes a description of the driver program and a description of the interface

to the lower level programs. It does not include a description of the detailed working of

the lower level routines themselves.

The drivers and example routines are provided in well-structured and indented Fortran 77.

The lower level routines were written in a private version of Ratfor [Kem75); the output

of the Ratfor processor is indented Fortran 77, but is not pretty. B2DE is single-precision.

1
Beginners to B2DE should first try one of the sample problems, using the interactive

driver. Its use is described in Section 3.2. To do a new problem, users must modify a main program, three subroutines, and a data file to describe the PDEs. They should read

Sections 3.1 - 3.3 and Section 3.4.1.

More experienced users may wish to change some of the USRxxx routines. They should

read all of Section 3.4.

Expert risers may wish to modify the drivers, the DRVxxx routines. They should read

Section 4.

Installation of B2DE on your computer and graphics terminals is covered in Section 5.

Section 6 contains several exzimple problems.

2. What B2DE Does

2.1 Mathematical Model

The current version of the B2DE solves systems of N PDEs, each of the form:

V- [ai{x,y,u,du/dx,du/dy)'7ui] = fi{x,y,u,duldx,du/dy) i 1,2,. ..,7V ( 1

)

in a region bounded by straight line segments and circular arc segments. On each segment

of the boimdary, each of the u* must obey either a nonlinear Dirichlet boundary condition

(type 1):

ffi(2T,y,u)=0 (2)

or a nonlinear normal derivative, of Neumann, boundary condition (type 2):

dtL

+ gi{x,y,u) = 0. (3)

On a given boundary segment, the type of boimdary condition may be different for different PDEs. This is a fairly general example of a system of elliptic PDEs arising from

conservation of a flux.

The numerical software was inspired by finite-element software of Bank amd Sherman

[Bank79, Bank82] and retains most of their philosophy.

At the heart of the package is a module which solves a system of linearized elliptic PDEs

using linear finite elements on a mesh of triangles. (Boundairy “triangles” may have one curved side, a circular arc segment.) The initial triangulation is am input to the package;

succeeding triangulations axe calculated adaptively by the package. For a given triangulation with M vertices, the calculation proceeds as follows.

The M finite-element basis functions {6,^} are linear on eawrh triangle; bm is 1 at vertex

m and 0 at all other vertices. For the tth PDE, the set of vertices on type 1 boimdaries is denoted by D{ (for Dirichlet). The solution is approximated by a sum of basis functions

M

ti.(a:,y) = X] “»m6m(x,y). (•*) m=1

2
The coefficients are determined by a Galerkin method [Stra73]; the error in the ith PDE

is made orthogonal to eax:h of the basis functions except those in J9»-:

jj

[- V • (a,- Vu,) + /,] = 0, m = 1, 2, . . . , M, m not in Di (5)

The remaining conditions on the a’s are that the type 1 boundary conditions hold exactly: = 0, m in Di. (6)

Equation (5) is integrated by parts to yield

jj

(a, Vu, • Vhm + fihm) + j

cagibm. = 0, (7)

where the single integral is over the type-2 parts of the boundary only. Equations (6) and

(7) comprise MN nonlinear equations in MN unknowns, the coefficients

The MN nonlinear equations are solved by an iterative process, a damped Newton’s

method. Let v be the Newton correction to u; t;,- has the same form as u»-. M

( 8) n»=l

To calculate v, replace u* by Ui + in (6) and (7) and linearize in v,-. The linearization

of (6) yields

N

ym)) ^ ^ •_ j-» fr\\ Qi J/m)) ^ r Pjm — Oj 771 in 2?,', (9)

and the linearization of (7) yields

jj{oiVvi ^bm+

N

Vu. . V6,„

y=i

dai dai dvi dai dvi

N r

y=i

duj ^ d{dujldx) dx d{dujldy) dx

}

dfi dfi dvj dfi dvj

duj^^ d{dujldx) dx d{dujjdx) dy

+

r ^ /‘E dgi dai dai dvj dai dvi

+ Si—«i + 9i -dv

--

; ^

+ g.- ' j=l ^ duj d{dujfdx) dx d{duj/dy) dy

3
( 10) Vu,- • V6„» + fibm Qi^m —

The right-hand side is the residual, r,m(u). Substituting for u and v from (4) and (8)

and doing the integrals numerically give a set oi MN linear equations for the correction

coefficients {/3im}- The matrix of the equations is the Jacobian matrix. For some problems,

the Jstcobian matrix may be singular, or nearly so. B2DE allows the user to “reguleLrize”

the Jacobian matrix by adding a constant to the “a” term in the Jacobizin matrix; this

is effectively like adding a multiple of the identity matrix to the Jacobian matrix. (The

residual is not changed, so that the solution foimd is the correct solution to the nonlinear

equations.)

B2DE allows the user to iterate on any combination of the PDEs, leaving the others fixed.

There is no guarantee of convergence if the initial guess is not good enough; how “good”

is “good enough” is problem dependent.

These linear equations are sparse; there are typically about 7N nonzero elements per row,

on average, and MN can be a few himdred or a few thousand. On the initial triangulation,

with MN a few hundred, a satisfactory way to solve these equations is directly, with

sparse Gaussian Elimination [Sher78|. On later triangulations, with MN up to a few

thouszmd, direct solutions take too much time and space; iterative methods are usually

better, especially multi-level iterative methods.

The linear equations on the coarsest level are solved by sparse Gaussian elimination. The

linear equations on the higher levels are solved by multi-level iteration [Bank79, Bank82,

Bran77]. The smoothing iterations provided by B2DE are Gauss-Seidel, minimum- residual

Gauss-Seidel, Block Gauss-Seidel using N equations per block (for the N unknowns at a

point), minimum-residual Block Gauss-Seidel, and Conjugate Gradients. The user specifies

r, the number of multi-level cycles, and m, the number of smoothing iterations in each

paxt of the cycle.

For systems of PDEs, the linear equation matrices are nonsymmetric, and may be difficult

to solve; there is no guarantee of convergence. Other smoothing methods would be desirable. Of the types in B2DE, Gauss-Seidel smoothing has been most satisfactory, even for

problems without theoretical reason to expect convergence.

Unless u is close to a solution of the nonlinear finite-element equations, u v may be

worse thzin u. A damped Newton’s method is used: replace u by u -t- A v, where A is chosen so that j] r(u + Av) ||<|| r(u) |1 where || r || is the Euclideatn norm of r. This avoids

divergence of the iterative process, but does not guarantee convergence.

In practice, the initial triangulation is usually too coarse to give the desired accuracy;

local refinement is needed to represent the solution accurately in part or all of the region.

After the package has converged to an approximate solution on the initial triangulation,

the approximate solution may be used to estimate the error and produce a new and finer

triangulation adaptively. There are various ways to estimate the error in each triangle.

For the present work, a local method is iised.

Consider a single triangle; for simplicity we suppose its vertices are 1, 2, and 3. The

4
approximate solution inside the triangle is

3 = X] (11)

m=l

We tentatively divide the triangle into four similar triangles by connecting the midpoints

of the three sides. On the four triangles, a potentially more accurate u,- is

= XI (12)

m=l

(The primes refer to the finer trizmgulation.) We keep the old values at the original 3

vertices; = aim,^ — 1»2,3. Values at the 3 new vertices axe foimd by solving the

linearized finite element equations for m = 4,5,6. The error in the triangle can be

estimated by the largest change at any of the points 4, 5, and 6.

This estimate is not accurate in absolute value, but the relative estimates in different

triangles are sufficiently accurate for refining the triangulation.

The refined triangulation is obtained by dividing the triangles with the largest estimated

errors, such &s all triangles with estimated error larger than 1/4 of the largest error, or the

worst 1/4 of the triangles. If the refined triangulation adds only a few more triangles, the

user may choose to “compress,” or merge the new triangulation with the previous one. For the refined triangiilation, the previously calculated u is probably a good guess at the

new u, so that few I'l’ewton iterations should be necessary. The hardest work is usually

done on the initial triangulation, the one with the fewest vertices.

2.2 Program Description

B2DE is written in a private version of Ratfor [Kem75]; the result is portable Fortran 77.

A few routines are unavoidably nonportable; these are described in Section 5. The PORT library [Fox78] machine constants are used to improve portability. All work

space is managed using the PORT library stack.

B2DE may be viewed as a three-layer program package. In the top layer are the programs

written by the user to describe a particular problem, and the USRxx programs which are user-modifiable. In the middle layer Bite programs callable by users, such as SOLxxx and

PLTxxx. In the bottom layer are programs which the user cannot see. The user does

not need to know anything about the bottom layer programs, and only needs to know

the calling interface to the middle layer programs. In particular, all the data structure of

B2DE is invisible to the user. B2DE comes with interactive driver programs, DRVxrx, which call middle- and top-layer

programs. The user is free to revise the driver programs.

5
3. Using B2DE with Supplied Drivers

3.1 Specifying Your Problem

To iise B2DE on your program, you need to write a main program and subroutines to

evaluate the coefficients o, /, and g and their first partial derivatives. The easiest way

to do this is to modify 2m existing program. See Appendix A and the sample programs

on the distribution tape for examples. The current DRVxxx driver assumes a data file to

specify the geometry of the region and a few parameters. See Appendix B and the sample

problems on the distribution tape.

Working space is managed inside B2DE using the PORT stack [Fox78]. Your main progrsm

must define the size of the stack and initialize it using iatkin, as in Appendix A. The size

of the stack determines the largest problem that can be solved.

Default versions of the USRxxx routines exist, but can be replaced by the user. See Section

3.4. For nonlinear problems, you will want to replace USRGSS to provide an initiai guess

at the solution.

Note that all of B2DE, including the DRVxxx routines, uses free-format reading. If you

make a mistake, such as putting a decimal point in an integer, or putting in a random

character, you will be informed and czm make a correction if you are running interactively.

If, for example, a program asks for four integers they can all be on one line, or on separate

lines. Blank lines are completely ignored. Input lines can have comments in them; a # or a

* begins a comment, which continues to the end of the line.

See Section 6 for the examples.

3.2 The Interactive Driver - Finding a Solution

DRVxxx programs zu*e drivers that make calls to subroutines in B2DE to solve the PDFs.

You may modify or replace any of these routines, but you should follow the model of the

default version carefully.

Your main program should be modeled on USRMAI. It will probably call DRVALL, which is the driver which hemdles the input phase joid then calls the solution driver. The default

version calls USRINl, then SOLINP, USRIN2, USRIN3, and finally DRVSOL. DRVALL assumes

that needed files have been opened and that the input/output values have been set by calls

to lOSET. Note that ausr, fusr, and gusr must be declared EXTERNAL; the names can be

anything you like, but ausr, fusr, amd gusr will be used in this manual.

subroutine DRVALL (ausr. fusr. gusr)

external ausr, fusr. gusr

DRVSOL is the driver for solving the PDEs. DRVSOL first calls TRMPLT to do any needed

initializing for plotting, then calls SOLREF to get started. The main part of DRVSOL is a loop in which an integer is read and one of nine actions is taken. Eight of the actions are calls to subroutines DRVPRT, DRVPLT, SOLITR, SOLREF, SOLCMP, SOLSAV, SOLD IF, and

SOLCHG. All th^e are described separately, and may be called by the iiser in a modified

driver. The ninth action is to quit; TRMPLT is then called to do any needed finishing off for

plotting.

subroutine DRVSOL (ausr, fusr, gusr)

6
external ausr, fusr, gusr

The remainder of this section is a description of the interactive driver, using a run from

Example 1. The responses from B2DE are those from print level 2. The interactive driver is menu-driven; user responses are one or more numbers. Except at

the top level, entering 0 as the answer to a question gets you to the next-higher menu.

On entering, B2DE prints the level (this will be 1 unless an old solution is being retrieved),

the time taien to set up the mesh (grid) , the number of vertices in the level 1 mesh (nv)

and the time taken to do initisd set-up (asmbnl).

Then the main menu is printed.

refsol entered at level 1 7 vertices

grid .033 sec. 16 vertices

asmbnl . 033 sec

l=print, 2*plot, 3=iter, 4=refine, 6=comp, 6=save, 7=diff, 8=end, 9=change

The possible responses will be described separately.

3.2.1 Printing

The DRVPRT driver program is called. Its main printing menu DRVPRT is l=8ummary, 2=timing, 3=storage, 4=vertices, 6=pde variables, 6=settings

Printing is to the terminal unless the 6=settings menu has been used to direct printing

and plotting to files. Each of the responses will be treated in order.

3.2.1.1 Summary

PRTSUM is czdled. For each level, the number of vertices and the number of triangles are

printed. For each PDE and each level for which refinement has been attempted, a measure

of the solution size, unorm, an estimate of the size of the error, enorm, and the estimated

number of correct significant digits in the solution are printed. The error and the estimated

number of digits aie only estimates; they are not guaranteed; they are almost certainly

not accurate.

n Ivl vertices triangles xmorm enorm digits 11 62 114 7.04E-01 1.04E-02 1.83

1 2 103 258

3.2.1.2. Timing

PRTTIM is called. The number of seconds taken in each of the major parts of B2DE

is printed, together with the percentage and the total, assembly is the time taken to

calculate Jacobian matrices, factor is the time taken to factor the level 1 matrices, mg

is the time taken to solve the linear equations to get the Newton step, rhs check is the

time taken to evaluate the residuzd to see if a Newton step is acceptable.

timing distribution

seconds percent

plot 0.000 0.00

7
assembly 2.200 41.51

factor 0.210 3.96

mg 0.180 3.40

rhs check 1.780 33.58

refine 0.930 17.55

total 5.300 100.00

3.2.1.3. Storage

PRTSTO is called. The maximum number of triangles and vertices are printed, together

with the number used so far. The maximum number depends on the size of the stack used.

103 vertices used.

2094 vertices allowed.

263 triangles used.

5583 tri2mgles allowed.

3.2.1.4. Vertices

PRTVER is called. The coordinates of each vertex are printed, together with the “parents”

of the vertex; the parents are the two vertices of the triangle edge which was bisected to

produce the vertex. (The original, user-defined vertices are their own parents.)

parents X y

1 1 1 -5.00000E-01 O.OOOOOE+00

2 2 2 O.OOOOOE+00 5.00000E-01

3 3 3 5.00000E-01 • 5.00000E-01

4 4 4 O.OOOOOE+00 O.OOOOOE+00

5 5 5 5.00000E-01 O.OOOOOE+00

6 6 6 O.OOOOOE+00 -5.00000E-01

7 7 7 -5.00000E-01 -5.00000E-01

8 2 4 O.OOOOOE+00 2.50000E-01

9 1 4 -2.50000E-01 O.OOOOOE+00

10 1 2 -3.53553E-01 3.53553E-01

11 3 4 2.50000E-01 2.50000E-01

3.2. 1.5. PDE variables

PDE values are printed. If a transformation is in effect, the transformed valueds are printed

when PDE 1 is asked for. PRTPDE is called. Its main menu is : l»on grid, 2=at vertices, 3=on grid as scaled integers

Values can be obtained on a regular rectangular grid or at the vertices. Grid values can be

printed as floating-point numbers or as scaled integers (xiseful for quickly getting an idea

of the solution). The three types will be discussed in order.

3.2.I.5.I. On Grid

A PDE is evaluated on a rectangular grid aligned with the i- and y-axes. If the magnification is 1, the grid is on a rectangle defined by the minimum and maximum i and y

vertex values. If the magnifications is greater than 1, the grid covers only a portion of this

rectangle.

8
The next menu is: enter pde no., l=value only, 2=gradient only, 3-both

Enter the PDE number, from 1 through npde, then 1, 2, or 3.

nx, ny

You are then asked for the number of grid points along x ajid along y. If the region

is not rectangular, or if its edg^ are not aligned with the axes, some grid points will

fall outside the region. These points have a large, machine-dependent number printed to

denote out-of-bounds; in these examples, the number is 1.7014lE-(-38.

A typical printing is:X

1 -5.00000E-01

2 -2.50000E-01

3 O.OOOOOE+00

4 2.50000E-01

5 5.00000E-01

6 -6.00000E-01

y

-5.00000E-01

-5.00000E-01

-5.00000E-01

-5.00000E-01

-5.00000E-01

-2.50000E-01

u

O.OOOOOE+00

O.OOOOOE+00

O.OOOOOE+00

1.70141E+38

1.70141E+38

4.03474E-01

3. 2. 1.5. 2. At Vertices

You axe asked which PDE to evaluate. Since the gradient is not continuous at the vertices,

it is not printed. •

enter pde no. A typical printing is:X

1 -5.00000E-01

2 O.OOOOOE+00

3 5.00000E-01

4 O.OOOOOE-^00

5 6.00000E-01

6 O.OOOOOE+00

y

O.OOOOOE+00

5.00000E-01

6.00000E-01

O.OOOOOE+00

O.OOOOOE+00

-5.00000E-01

u

l.OOOOOE+00

l.OOOOOE+00

6.68726E-01

7.04038E-01

O.OOOOOE+00

O.OOOOOE+00

3. 2. 1.5.3. On Grid as Sc^ed Integers

A PDE is evaluated on a rectangular grid aligned with the x- and y-axes. If the magnification is 1, the grid is on a rectangle defined by the Tninirmim and maximum x and y

vertex vedues. If the magnifications is greater than 1, the grid covers only a portion of this

rectangle.

The next menu is: enter pde no., l=value only, 2“gradient only, 3= both

Enter the PDE number, from 1 through npde, then 1, 2, or 3. nx, ny

9
You axe then asked for the number of grid points along x and along y. If the region is not rectangular, or if its edges axe not aligned with the axes, some grid points will fall outside the region. These points have 9999 printed to denote out-of-bounds. Other points

are scaled to lie from -100 to 100. A typical printing is: X limits -1.66667E-01 1.66667E-01

y limits -6.00000E-01 -1.66667E-01

pde 1 scaled values from -6.085E-08 to 5.492E-01 scale 5.492E-01

90 90 91 93 95 97 100

78 78 80 82 85 89 92

64 65 67 70 75 80 85

49 50 54 58 65 71 77

34 35 39 44 54 62 70

17 19 21 28 42 55 65

0 0 0 0 9999 9999 9999

3.2.1.6. Settings

The next menu is l=color, 2=erase. 3=magnify, 4=window, 6*file print /plot, e^transform

The same menu is available from the plotting menu. All six menu choices afiFect plotting,

but only choices 3, 5, and 6 affect printing. This menu is discussed in detail in Section

3.3.8, under plotting.

3.2.2 Plotting

Plotting is covered in great detail in Section 3.3.

3.2.3 Iterating

Using the current mesh, do one or more Newton iterations on the nonlinear finite element

equations which represent the PDEs. In each iteration, a Jacobism matrix emd a righthzmd-side are calculated; the solution of these equations is the Newton correction. If the

mesh level is 1, the linear equations are solved by sparse Gaussian elimination. If the mesh

level is greater than 1, the linear equations are solved by a multi-level iteration, with the

level 1 equations solved by spjirse Gaussian elimination.

If the PDEs are linear, one iteration is suflBcient if the linear equations are solved sufficiently

accurately. The prompt message is enter niter, ieps, relerr. abserr. jsmax. param

niter is an integer, the maximum number of Newton iterations to be done. Fewer

will be done if the accuracy requirements are met.

ieps is an integer which governs the accuracy with which the linear equations arc solved. Let || r H be the Euclidean norm of the initial right-hand-side; this is also

the residual of the nonlinear equations and the initial residual of the linear equations.

The multi-level iterations are stopped when the new residual is

|| r j| or less, or when maxr multi-level iterations have been done, maxr is set by subroutine USRIN3.

10
Tlie value of ieps is immaterial if the level is 1, since level 1 equations axe solved by

sparse Gaussian elimination.

If ieps is too small, the linear equations will not be solved accurately enough, and more

iterations may be needed. If ieps is too large, more work may be used in solving the

linear equations than is justified by the accuracy of the approximate solution, relerr

and abserr are real numbers which govern the accuracy with which the nonlinear

equations are solved. Newton iterations are stopped when the residual is small enough,

when II r 1|<= (relerr || u || -fabserr), where || u || is the Euclidean norm of the

current approximate solution.

j smax is an integer which governs the damped Newton solution. If v is the Newton

correction, the new approximate solution is u + tv. First t = to is tried, where to is the value returned by USRSTP, then to/2, then to/4, and so on, until the new residual

is small enough; that is, imtil || r(u-f-tv) |1<|| r(u) ||. A total of jsmax tries are made.

param is a real parameter for the convenience of the user; it is put into a common

block and can be accessed by the user. It may be useful in, for example, USRREG. The

common block is COMMON /USRCOM/ PARAM

If the nonlinear iterations have not converged by the niter-th iteration, DRVSEE is called

to give the user an opportimity to look at the Newton correction before the correction has

been trated. The menu is l=print, 2=plot, 3=continue

Menu choices 1 and 2 are exactly the same as choices 1 amd 2 of the main menu, except that

the user prints or plots the current Newton correction rather than the current approximate

solution. The current magnification, window, and color are used. Use menu choice 3 to go

on to test the Newton correction.

B2DE has provision for user-supplied adjustments. For example, the PDEs may depend

on the solution to some nonlinear equation formally separate from the PDEs. After u has

been replaced by u -t- tv, but before the residual is calculated, USRADl is called. After a

satisfactory value of t has been found, and after the residual is calculated, USRAD2 is called.

3.2.4 Refining

Using the current approximate solution, estimate the error in the solution and adaptively

refine the mesh. B2DE first displays the current level number, the number of vertices, and

a crude guess at the maximum error in the solution for each PDE. A typical example is: refsol re-entered at level 2 26 vertices

asmbnl . 050 sec

1 estimated maximum error 8.01E-03

B2DE then gives you an opportimity to plot the estimated error for each PDE. When

you are finished, if npde is 2 or more, B2DE gives you an opportunity to plot the total

estimated error for all the PDEs combined.

B2DE then shows the estimated Euclidean error for each PDE, shows how many triangles

11
will be refined because the local estimated error is high enough, and then refines. The

estimated error and the number of correct digits axe almost certainly wrong, but are indicative of the progress of the solution from level to level.

If any of the weights returned by USRWTS is negative, B2DE asks for the weights to use and asks if you want to refine.

level 1 pde 1 unorm, enorm 7.06E-01 2.86E-02 eat. digits 1.39

min error / mzoc error 1.44 percent

4 of 30 triangles ( 13.33) above 60.00 percent

9 of 30 triangles ( 30.00) above 25.00 percent

8 of 30 triangles ( 26.67) above 37.50 percent

B2DE tells you the number of vertices at the new level.

refine .700 sec. level 2 52 vertices

3.2.5 Compressing

The vertices of the highest-level mesh are merged with those of the next-to-highest- level

mesh, and the level number is reduced by 1. This option is useful in cases where a refinement generates only a few new vertices, as may

happen in problems with a singularity in the solution.

3.2.6 Saving the solution

The current approximate solution is saved on a file. The file must have been opened by

the user before calling B2DE. Let n = ioget(8).

If n > 0, the information is written as unformatted (binary) information on logical file \mit

n. If n< 0, the information is saved as formatted information on logical file unit -n. The

prompt message is enter the number of significant figures to save

After the user enters the number of significant figures to save, B2DE constructs an appropriate format and writes the file. Formatted files tadte more space to write, but can be used to transport data from one computer to another.

3.2.7 Looking at the Difference

Look at the change in the current solution since the last time the mesh was refined, or

since the start of the current session with B2DE, whichever is later. DRVSEE is called.

The menu is l=print , 2=plot , 3*continue

Menu choices 1 and 2 are exactly the same as choices 1 aind 2 of the main menu, except

that the user can print or plot the ciirrent difference rather than the ctirrent approximate

solution. Use menu choice 3 to go back to the main menu.

3.2.8 Quitting

12
Quit.

3.2.9 Changing Parameters

Change user-defined parameters, internal B2DE parameters, or change which PDE variables are being iterated on.

First USRCHG is called; in USRCHG, the user may change any xiser-defined parameters as

desired.

If npde is 2 or more, the user may keep the values of one or more of the PDE variables

fixed, and calculate Newton corrections for the remaining PDE variables. If npde is 2 or more, B2DE next prints

PDE variables - 0 for fixed, 1 for varying

2nd prints the current status of each of the npde PDE variables,

enter new values

The user must enter the new values as integers: 1 means the PDE is to vary, and einy other

integer means the PDE is held fixed. At least one 1 must be entered.

The next prompt is enter 1 to change pareimeters , 0 to return

If 1 is entered, the next prompt is parameters: maxm, maxr, ittype, ithrsh, iprsw

and the current values are printed out. Enter the new values, all integers. Refer to Section

3.4.9, USRIN3, for the meanings of the parameters. m2Lxm, maxr, and ittype all govern

multi-level iteration, ithrsh governs refinement, iprsw governs the amoimt of printing as the solution progresses.

3.3 The Interactive Driver - Graphics

Seeing pictures is necessary if the user is (l) to understand what is happening during the

iterative process leading to a solution and (2) to understand the properties of the solution.

B2DE gives the user ample opportunity to msike plots. There are standard plots available

of several kinds. The variables plotted axe either a PDE variable or the magnitude of its

gradient. Since the author cannot guess everything the user might want to plot, B2DE

allows the user to define problem-dependent plot variables (Section 3.4.12). Plots may

be made to the user’s graphics terminal, or the vectors may be written to a file for later

processing.

Most of the graphics progrsim code is independent of the computer used and of the graphics

output device. The dependent parts are treated in the installation part of the manual.

Section 5.

Five kinds of plots are available: triangles, contours, surfaces, profiles, and flow lines. Each

plot type is done by a subroutine which may be called directly from a xiser’s program. Each

plot type also has an interactive driver which may be called directly from a user’s program.

There is &n overall interactive plot driver.

13
The interactive driver is DRVPLT. The user may wish to modify or replace this progrsim.

DRVPLT asks interactively for user input, then calls any or all of DRVTRI, DRVCON, DRVSUR, DRVPRO, DRVFLW, and DRVSET. The user may wish to modify or replace 2iny of these

programs.

Most of the DRVxxx programs call the corresponding PLTxxx programs: PLTTRI, PLTCON,

PLTSUR, PLTPRO, and PLTFLW. The calling sequences for the PLTxxx programs are described

in Section 4.3; PLTxxx programs may be called directly by the user, but are not meant

to be changed by the user. Some PLTxxx programs require work space, which is obtained

from the stack; if insufficient stack space is available, the plot is skipped.

DRVSET calls any of six sub-driver programs: DRVCOL, DRVERA, DRYMAG, DRVWDW, DRVFPL, and

DRVSHF. The DRVxxx programs call the corresponding SETxxx programs: SETERA, SETCOL,

SETMAG, SETWDW, SETFPL, and SETSHF. The calling sequences for the SETxxx programs are described in Section 4.1; SETxxx progrzuns may be called directly by the user, but are not

meant to be changed by the user.

Plots can be drawn to the full screen or plotting surface, or can be drawn in a subwindow.

Contour plots, flow line plots, and triangle plots can be magnified.

All plots axe of vectors only, with no plot labeling. If the hardware \ised supports it, erasing

between plots is optional. If the hardware used supports it, the line type or color of the

plots can be chamged.

The following description of plotting assumes the overall interactive driver, DRVPLT. All

the DRVxxx programs are menu-driven and us*e a free-format input package. As in the rest of B2DE’s interactive drivers, responses of 0 return the user to the next higher-up menu.

Erroneous responses are dealt with reasonably; error messages should be self-explanatory.

The top-level menu of DRVPLT is: l^triangles, 2=contour3, 3=surface, 4=profilea, 6=flow lines, 6=3ettings

The response should be a non-negative integer. A response of 0 meems to go to the next

higher menu. Other responses outside the indicated range are ignored. The responses will

be treated in order.

3.3.1 Triangle plots

The DRVTRI driver is called. The next menu is : enter level

Any existing level of the triangulation can be plotted. The current magnification and

window are used. The plot is in the current color.

After the plot is done, the program may wait for the user’s signal. The signal, if any, is that required by the terminal-dependent program, PLTUTL. After the signal, the time taken

is displayed, and the

enter level

menu comes up again. To get back to the previous menu, enter

0 0

14
3.3.2 Contour Plots

The DRVCON driver is called. The next menu is : pde no., l=f;mction/2=ab8(grad) , l=fill/2=no fill

Enter the PDE number, from 1 through iV, the number of PDEs in the problem. Enter 1

for contours of the PDE value, 2 for the absolute value of its gradient. Enter 1 to fill in

the contours with different colors in the current palette, or 2 to just draw the contours in

the current line color. The current magnification and window are used. If fill is chosen,

the next menu is : start with color number

Enter the first color of the color palette, which may be from 0 on up. The space between

the bottom two contours will be filled with this color; succeeding spaces will be filled with

successive colors. The next menu is : number of contours

The number of contours is limited to 101 by PARAMETER (MAXCON = 101) in program

DRVCON. The magnification is treated as for triangle plots.

If nc, the number of contours, is positive, the program finds the maximum and minimum

values of the plot variable and spaces the contours equally between, but not including, the

extreme values. Both the extreme values and the contour values are displayed.

If nc is 0, the extreme values are displayed, but no plot is done.

If nc is negative, the user is prompted for extreme values for contours. There are -nc

contours, equally spaced between and including the extreme values.

If no fill is chosen, the next prompt is : number of contours

and the replies are as above.

After the plot is done, the time taken is displayed, and the

pde no., l»function/2»abs(grad) , l=fill/2=no fill

menu comes up again. To get back to the previoiis menu, enter

0 0 0

3.3.3 Surface Plots

The DRVSUR driver is called. The current color and window are used. The magnification

is 1. The next menu is : pde no., l=function/2=abs(grad)

Enter the PDE number, from 1 through N. Enter 1 for a plot of the PDE value, 2 for a

plot of the absolute value of its gradient. The next menu is : level, l=surface triemgles, 2=surface contours, 3=surface grids

Enter the triangulation level desired for the plot, and the type of stirface plot desired. The

absolute value of the level is used. In all three types, the function is plotted as a parallel

15
projection plot (no perspective); the height of the surface above or below the x-y plane

represents the function value. The surface can be viewed from any angle. Hidden lines are omitted unless the level is negative; then hidden lines are plotted. Plotting without

removing hidden lines requires much less computation, and can be much faster.

The three types will be discussed in order.

3.3.3.1 Surface triangles

The surface is defined by the triangles of the specified level. The next menu is : enter viewing direction — x, y, z The plot is displayed as if seen by a viewer looking towards the surface along a line

extending from the origin through the point (x, y, 2). After the plot is done, the time

taken is displayed, and the

pde no., l=fxmction/2=abs (grad)

menu comes up again. To get back to the previous menu, enter

0 0

3.3.3.2 Surface contours

The surface is defined by the contour lines drawn on triangles of the specified level. Triangle

boundaries are not shown. The first menu is again : pde no., l=function/2=ab8(grad)

and the responses are as before. The next menu is : number of conto\irs

Responses to this menu are like those for the planar contour plot, except no magnification is allowed. The number of contours is limited to 202 by PARAMETER (MAX » 202) in program

DRVSUR.

The next menu is : enter viewing direction — x, y, z Responses are as before. After the plot is done, the time tadcen is displayed, and the

pde no., l=function/2=abs(grad)

menu comes up again. To get back to the previous menu, enter

0 0

3.3.3.3 Surface grids

The surface is defined by grid lines parallel to the x zind y aoces, drawn on triangles of the

specified level. Triangle boimdaries Me not shown. The first menu is again

pde no., l=function/2=abs(grad)

and the responses are as before. The next menu is number of x lines, n\imber of y lines

16
The X lines £ire equally spax:ed along the x-axis, from the minimumx-value to the maximum

X-value; the y lines are similarly spaced along the y-axis. The total number of grid lines, x

plus y, is limited to 202 by PARAMETER (MAX = 202) in program DRVSUR. The next menu

is: enter viewing direction — x, y, z Responses are as before. After the plot is done, the time taken is displayed, and the

pde no., l=function/2=abs(grad)

menu comes up again. To get ba^k to the previous menu, enter

0 0

3.3.4 Profile Plots

The DRVPRO driver is called. The current color zmd window are \ised. The next menu is : enter pde no. , number of profiles

Enter the PDE number, from 1 through iV, and the number of profiles desired. A profile

is a slice through a surface plot, viewed from the side. The slices Eire parallel and equally

spaced. The next menu is : enter indices of profile plane (s) x, y

The slices are parallel to a line running from the origin to the point (x, y). Viewing is from the direction as defined for surface plots, (y, -x, 0), perpendicular to the slices. After

the plot is done, the time taken is displayed, and the

enter pde no., number of profiles

menu comes up again. To get back to the previous menu, enter

0 0

3.3.5 Flow Line Plots

The DRVFLW driver is called. Flow line plots attempt to describe the fiux, or vector field, a, Vui. The user specifies a line segment, or set of connected segments; the positive and

negative fiuxes across the line are calculated; points along the line are spaced evenly in

the fiux; and flow lines 2U“e drawn both ways from these points. Flow lines axe terminated

at boimdaries of the region or when the fiow lines would re-enter a triangle. The current

magnification and window are used. The plot is in the current color.

Since the solutions are composed of linear finite elements on triangles, Vu» is constant

within each triangle. Gradients tend not to go smoothly to zero when they should, but

overshoot instead; this makes the fiow lines stop.

The next menu is : enter pde no., number of points, number of flow lines

The PDE number can be from 1 to N. The number of points is the number in the set of

connected line segments, two or more. The next menu is : enter # x-values followed by # y-values

17
where # is the number of flow lines just entered. The total flow in a positive and in a negative directions across the broken line is printed, and then the flow lines are drawn.

After the plot b done, the time taken is displayed, and the

enter pde no., number of points, number of flow lines

menu comes up again. To get back to the previous menu, enter

0 0 0

3.3.6 Settings

The DRVSET driver is called. Several values can be set; these apply to future plots until

the settings are changed. The menu is l=color. 2=erase. 3=magnify. 4=window. 5=file print/plot, 6=transform

Setting the Color

The DRVCOL driver is called. Some terminals can support vectors drawn in different colors or line types. The menu is : l^change color number. 2=:change color palette

If 1 is entered, the next menu is: enter new color number

SETCOL is called, which calls TRMCOL, and then the settings menu appears.

If 2 is entered in response to

l*chaoige color number, 2=change color palette

the system-dependent subroutine TRMPAL is called. Then the settings menu appears.

Setting the Erase Flag and Erasing the Screen

The DRVERA driver is called. Some terminals support optionad erasing between plots. If the erase flag is ON, the screen is erased between plots; if the erase flag is OFF, the screen

is not erased between plots. You might want to plot, for instance, s\irface contours on the

same plot as a surface grid, or planar contours of more than one PDE. The menu is: l=set erase flag off. 2=set erase flag on. 3=erase screen

If option 3 is chosen, the screen is erased immediately, but the erase flag is not changed.

There are no lower menus. After the desired option is chosen, SETERA is called and the

settings menu appears.

Setting the Magnification

The DRVMAG driver is called. B2DE can make magnified triangle, flow line, and contour

plots. The parameters necesary are the magnification and the point which will be the

center of the magnified plot. The magnification also affects printing of PDE variables on

a grid, either as real values or as scaled integers. The initial menu is: enter magnification

If the magnification entered is 1 or less, the magnification is set to 1 and you are returned

to the main plot driver menu. If the magnification is greater than 1, the next menu is:

18
enter x-center am.d y-center (0 to 100)

An unmagnified plot fits in a rectangle aligned with the x- and y-axes, defined by the

TniniTTmm £ind maximum x and y vertex values. For magnifications greater than 1, the

center of the magnified plot is defined with respect to this rectangle; the center is defined by

the percentage distance from the minimum to the maximum coordinates of the rectangle.

For plotting, the center may be adjusted to fill as much of the screen as possible. For

example, if the center is given as 100 100, the upper right comer of the rectangle will be

plottd to the upper right comer of the screen rather than to the center of the screen. SETMAG is called and the settings menu appears.

Setting the Window

The DRVWDW driver is called. B2DE can draw its plots in the entire plotting surface as defined by TRMSCL, or can use a subwindow. The prompt is enter x-left, x-right, y-bottom, y-top (from 0 to 100)

The numbers entered are interpreted aa percentages of the full window. For exaonple,

0 50 50 100 means the upper-left quarter of the full window. SETWDW is called and the

settings menu appears.

Printing and Plotting to Files

The DRVFPL driver is called. Printing and plotting can be done to the screen of the terminal

or to a file. One fiag governs both printing and plotting. The prompt is l=plot/print to screen, 2=plot/print to file

SETFPL is called and the settings menu appears.

Users may write their own programs to read the plot file and produce plots on other

devices. The format of the file is as follows. The first line is produced by the terminaldependent program PLTUTL, and may be changed by the user. The default version writes a line containing the minimum and maximum integer values produced by the plot,

xmin ymin xmax ymax

in (417) format. These values come from the terminal-dependent program TRMSCL, eind

may be changed by the user; the default values are 500, 9500, 500, and 7250.

The plot vectors are next in the file. Each is on a separate line,

xl yl x2 y2

in (417) format. The end of the plot is marked by a plot vector

-1 -1 “1 -1

which is produced by TRMPLT, and which may be changed by the user. Transforming a Solution

The DRVSHF driver is called. The user may specify transformations of PDEs and print or

plot the transformed values. The prompt is enter index of transformation

19
SETSHF is called, which calls USRSHF, and the settings menu appears.

3.4 Easy Modifications - USlbooc Programs

Default routine USRxxx aire supplied; these routines are designed to be modifiable by the

user. They provide ezisy ways to alter the performance of B2DE without going inside

B2DE. Users should read the default versions carefully before modifying the programs.

3.4.1 USRMAI is the Fortran “main” program. USRMAI declares storage, sets up for input

and output (I/O), and calls a driver routine. The a, /, and g programs should be declared

EXTERNAL in this program. The stack must be declared and initialized in the main program.

1/ O values must be set by calling lOSET for all nine values. These values may be retrieved

by calling lOGET; B2DE routines needing to do I/O get file numbers in this way.

io(l) is the file number to be used for reading from the terminal.

io(2) is the file number to be used for writing to the terminal.

io(3) is the terminal type; it is mainly used by graphics programs. If io(3) is nonpositive, the session is assumed to be batch, not interactive. All plots will be made

to files, and FREEIN, the free-format input routine, will abort if am input line is bad.

io(4) is the file number to be used for reading input data; it is only used in USRxxx

programs.

io(5) is the file number to be used for writing plot data, if any platting is done to

files rather than to the terminal.

io(8) is the file number to be used for writing print data, if any printing is done to

files rather than to the terminal.

io(7) is the file number to be used for reading the old solution data, if an old solution

is to be read. If io(7) is positive, the solution data is assumed to have been written

as unformatted data. If io(7) is negative, the solution data is assumed to have been

written as formatted data.

io(8) is the file number to be used for writing the new solution data, if a solution

is to be saved. If io(8) is positive, the solution data will be written as unformatted

data. If io(8) is negative, the solution data will be written as formatted data.

io(9) is the file number to be used by B2DE for writing scratch data. Data is written

to and read from io(9) as unformatted data.

Data files corresponding to file numbers io(4), io(5), io(6), io(7), and io(8) should

be opened if there is any chance they will be used. File io(9) is always used and should

be opened.

DRVALL or another driver should be called.

Any files opened should be closed.

3.4.2 USRADl is called to make ainy needed adjustments after a tentative Newton step has

been taken, but before the new residual is calculated. The default version does nothing.

subroutine USRADl

20
3.4.3 USRAD2 is called to make any needed adjustments after a Newton step has been

accepted, after the new residual is calculated. The default version does nothing.

subroutine USRAD2

3.4.4 USRCHG is called when the 9=change selection from the main menu is taken. It should return .TRUE, if the changes made affect the residuals of the PDEs so they must

be recalculated; it should return .FALSE, if the residuals need not be recalculated. The

default version returns . FALSE

logical function USRCHG (junk)

junk is a dummy parameter.

3.4.5 USRGET is called when an old solution is being read in, to read anything that was

saved by USRSAV. The default version does nothing. USRGET should read values in the

same order and with the same formats as USRSAV writes. USRGET is always called with

isin equal to io(7).

subroutine USRGET (isin)

integer isin

isin is an input, the hie number to read from.

3.4.6 USRGSS is called to get the initial guess at the solution. The default version guesses

0; this is adequate for linear problems, but is probably not adequate for most nonlinear

problems. USRGSS Is C2klled for each point, in the initial mesh.

subroutine USRGSS (npde, x, y, i. u)

integer npde. i

real x, y, u(npde)

npde is the number of PDEs.

X and y are the coordinates of the point.

i is the number of the user-supplied triangle in which the point lies,

u is an output, the guess at the solution at the point.

3.4.7 USRINl is the first input routine. It returns six values. The default version calls

I0G£T(4) to get a file number, and calls FREEIN to read the values from the file. USRINl

is called by DRVALL; if DRVALL is replaced by the user, USRINl may not be needed.

subroutine USRINl (oldsol, nt. nv. nc, nb. npde)

logical oldsol

integer nt , nv . nc , nb . npde

oldsol should be .TRUE, if an old solution is to be read in, and .FALSE, if a new

solution is to be started from the beginning. If oldsol is . TRUE . , the remaining five

values do not need to be returned; the old values sure used by B2DE.

nt is the number of triangles in the user-supplied coarsest mesh.

nv is the number of vertices in the user-supplied coarsest mesh.

nc is the number of curved boundary edges in the user-supplied coarsest mesh.

21
nb is the number of boundary edges in the user-supplied coarsest mesh,

npde is the number of PDEs.

3.4.8 USRIN2 is the second input routine. It is called to obtain the geometry of the user- supplied coarsest mesh. It is called only if oldsol is . FALSE . . The default version calls

lOGET (4) to get a file number, and calls FREEIN to read the values from the file.

subroutine USRIN2(nt, nv, nc, nb, npde, Iref, vxin, vyin,

* xmid, ymid, itnode, ibdry)

integer nt, nv, nc, nb, npde

real vxin(nv) , vyin(nv) , xmid(nc) , ymid(nc)

integer ibdryCnb, npde+4) , itnodeO, nt) , Iref(nv)

nt, nv, nc, nb, and npde are all inputs; they are as described under USRINl. The

remainder of the arguments are outputs.

Iref (j) is the initial refinement level for vertex j.

vxin(j) is the x-coordinate of vertex j.

vyin(j) is the y-coordinate for vertex j.

xmid(j) is the x-coordinate for the midpoint of curved edge j, if nc > 0. ymid(j) is the y-coordinate for the midpoint of curved edge j , if nc >0.

itnode (l,k), itnode (2, k), and itnode (3, k) are the three vertices for triangle k.

ibdry (m, n) contains data for boimdary edges m. ibdryCm, 1) must be m. ibdry (m, 2) is the number of the tri£Uigle containing the boimdary edge, in the range

1 to nt.

ibdry (m, 3) is the edge number of the boimdary edge within the triangle. It must

be 1, 2, or 3; the first edge is opposite the first-named vertex in the itnode array for

the triangle.

ibdry (m, 4) is 0 if the edge is straight, and is the curved edge number, in the range

1 to nc, if the edge is curved.

ibdryCm, n+4) is the boundary condition type for the n-th PDE on this boundary

edge; it should be 1 for Dirichlet-type boundary conditions amd 2 for Neumann-type

boundairy conditions.

3.4.9 USRIN3 is the third input routine. It is called to obtain parameters for the solution

process. It is called only if oldsol is .FALSE. The default version calls I0GET(4) to get a file number, and calls FREEIN to read the values from the file. All the airguments are outputs.

subroutine USRIN3(mauan, maxr, ittype, Ivlref , Ivlusr, numlvl,

* ithrsh, iprsw, iwork)

naxm is the “m” parameter for multi-level iterations [BaLnk79|.

22
mayr is the “r” parameter for multi-levei iterations [BankTO].

ittype is the type of smoothing iteration.

0: point Gauss-Seidel.

1: conjugate gradient.

2: block Gauss-Seidel.

Ivlref is the maximum pre-refinement level assigned to any vertex.

Ivlusr is the maximiim pre-refinement level used in USRTST.

numlvl is the number of mesh levels allowed,

ithrsh is the threshold percentage for mesh adapting,

ithrsh = 0 : do uniform refinement.

ithrsh > 0 : refine triangles whose estimated error is at least this percent of the

largest estimated error. ithrsh < 0 : refine triangles whose estimated error puts them in the worst

abs (ithrsh) percent of all triangles.

iprsw determines how much printing is done as the solution is being found.

0 : as little as possible.

1 : more.

2 : reasonable for interactive use.

3 : far too much.

iwork is the number of words of stack spew:e to use for solving. B2DE will use the

smaJlest of iwork, the value initialized by USRMAI, and the actual available space.

3.4.10 USRREG should return the regularization parameter for the n-th PDE. In calculating

the Jacobian matrix, B2DE uses a„ + USRREG (n) to allow for regularizing nezirly-singular

Jacobian matrices. The default version returns zero. real function USRREG (n)

integer n

common / usrcom / param

real param

3.4.11 USRSAV is called when a solution is being saved, for the user to save necessary

problem paurameters not known to B2DE. The default version does nothing. USRGET should

read values in the same order and with the same formats as USRSAV writes. USRSAV is always

called with isout equal to io(8).

subroutine USRSAV(isout)

integer isout

isout is the file number to write to. 3.4.12 USRSHF is called to “shift,” or transform PDE variables for plotting. All the argu23
ments are inputs except unew, which is an output. The default version does nothing.

subroutine USRSHFCixtype , npde, nv, u, unew, vx, vy)

integer ixtype, npde, nv

real u(npde,nv), unew(nv), vx(nv) , vy(nv)

ixtype is the index of the transform type.

npde is the number of PDEs.

nv is the number of vertices in the current finest mesh.

u is the solution on the finest mesh. vx and vy are the coordinates of the vertices.

unew is the transformed solution. It may be any desired function of the PDE variables,

the coordinates, or anything else.

3.4.13 USRSTP provides the user an opportunity to avoid problems such as overfiow or the

square root of a negative number, by restricting the maximum fraction of a Newton step

that can be used. B2DE will try to evaluate residuals at (uold -f- t uatep), where t is the

value returned by USRSTP. USRSTP should return a number greater than 0 and less than or equal to 1. All the arguments are inputs. The default version returns 1. real function USRSTP (npde, nv, vx, vy, uold, ustep)

integer npde. nv

real vx(nv) , vy(nv) , uold (npde, nv) , ustep(npde, nv)

npde is the number of PDEs.

nv is the number of vertices in the current finest mesh, vx and vy are the coordinates of the vertices,

uold is the ciirrent solution,

ustep is the tentative Newton step.

3.4.14 USRTST allows pre-refinement of the initial mesh in addition to that provided by

the array Iref in USRIN2. After the pre-refinement scheduled by Iref , USRTST is called

for each triangle, and should return .TRUE, if the triangle is to be refined at least one more time. All triangles are tested until no more triangles are refined. All the arguments

are inputs. The default version always returns . FALSE

logical function USRTST (x, y, Ivltri, level 1)

real x(3) , y(3)

integer Ivltri , levell

X and y are the coordinates of vertices of the triangle.

Ivltri is the level of refinement of this triangle so far.

levell is the maximum level of pre-refinement allowed, as specified by USRIN3.

3.4.15 USRWTS is called when the mesh is being refined; it should return the relative weight

for the n-th PDE. If any of the weights is negative euid the run is interewrtive, the driver

asks what the weights should be, ignoring the values returned by USRWTS. The default

24
version always returns 1. real function USRWTS(n)

integer n

n is the PDE number.

4. Writing Your Own Driver for B2DE

The routines callable by the user axe described in several groups: DRVxxx, SETxxx, SOLxxx,

PLTxxx, PRTxxx, TRMxxx, USRxxx, and miscellaneously-named routines.

DRVxxx routines are driver programs, and can be modified or replaced by the user. DRVALL

and DRVSOL are described in Section 3.2, and DRVPLT is described in Section 3.3.

DRVSEE can be used to inspect the solution

subroutine DRVSEE (ausr)

external ausr

DRVSEE provides the menu discussed in Section 3.3.3:

Imprint , 2=plot , 3=continue

SOLxxx routines do the solution of the PDEs. They are described in Section 4.2.

PLTxxx routines do plotting. They are described in Section 4.3.

PRTxxx routines do printing. They are described in Section 3.2.1.

SETxxx routines set parameters used for plotting. They are described in Section 4.1.

TRMxxx routines are terminal-dependent programs having to do with plotting. They axe described in Section 5.2.

USRxxx routines provide the user with opportunities to do things without requiring detailed

knowledge of the interior of B2DE. They axe described in Section 3.4.

4.1 SETxxx Routines

4.1.1 SETCOL sets the color or line type for plotting. It saves the input value, then calls

the terminal-dependent routine TRMCOL; TRMCOL is described in Section 5.2.1.

subroutine SETCOL (newcol)

integer newcol

newcol is the new color or line type.

4.1.2 SETERA either erases the screen or sets the fiag that determines whether the screen is erased before a new plot is drawn. On some terminals, it may not be possible or desirable

to plot without erasing the screen. subroutine seteraCif lag)

integer iflag

If iflag is 1, the screen will not be erased before each plot.

If iflag is 2, the screen will be erased before each plot.

If iflag is 3, the screen is erased immediately.

25
Other values for if lag do nothing.

4.1.3 SETMAG sets the magnification emd the x- and y-centers for the magnified plots.

Magnification applies to triangle plots, contour plots, and fiow line plots. It also applies

to printing of PDE values on a grid.

subroutine setmagCxymag, xctr, yctr)

real xymag, xctr, yctr

xymag is the magnification. If xymag is less than 1, B2DE uses a magnification of 1. xctr and yctr are the coordinates of the center of the image to be magnified, in

percentages of the full window size. For example, if xctr and yctr axe each 50, the

center of the image will be magnified.

4.1.4 SETWDW sets the window in which subsequent plots will be drawn.

subroutine setwdwCxmin, xmax, ymin, ymax)

real xmin, xmax. ymin. ymax

xmin and xmax aure the x-coordinates of the window, in percentages of the maxim\im

window size as given by TRMSCL. Similarly, ymin and ymax axe the y-coordinates of

the window. For exatmple, if xmin amd ymin are eau:h 50, and xmax aind ymax are each

100, the plot window will be the upper right quarter of the plotting area.

4.1. SETFPL sets a fiag for plotting and printing.

subroutine SETFPL (flag)

logical flag

If flag is .TRUE., subsequent printing and plotting will be done to a file.

If flag is .FALSE., subsequent printing and plotting will be done to the terminal.

4.1. SETSHF is used to “shift”, or transform, the first PDE variable for inspecting a user-defined function of the PDE variables, and to shift it back.

logical function SETSHF (if lag)

integer if lag

If if lag is not 0, SETSHF calls USRSHF to get new values to put in place of the values

of the first PDE, exchzmges the values, and returns for plotting or printing, if lag is passed to USRSHF for use as an index.

If if lag is 0, SETSHF restores the values saved for the first PDE.

4.2 SOLxxx Routines

SOLxxx routines are called by driver programs to solve the PDEs. They are described

individually below.

4.2.1 SOLCHG is used to change parameters during a solution. There are three types of

parameters. First, USRCHG is called for changing user-supplied pexameters. Second, if there are two or more PDEs in the problem, one or more PDE vEiriables can be held fixed during

the nonlinear iterations. Third, B2DE parixneters maxm, maxr, ittype, ithrsh, and iprsw

can be changed.

26
subroutine SOLCHG

4.2.2 SOLCMP is used to compress, or merge, the top two levels of mesh into a single mesh.

Lower level meshes are not affected. SOLCMP is useful if a refinement did not add many

new vertices to the mesh. There are no inputs or outputs.

subroutine SOLCMP

4.2.3 SOLDIF is used to inspect the “difference” in a solution; this is the difference between

the starting solution values (initial guess, read-in values from old solution, or interpolated

solution from next coarser mesh) on the current finest mesh ajid the current solution

values. Norms of the changes are printed, and DRVSEE is called for possible printing or

plotting. SOLDIF is useful for inspecting the changes in the solution. The only input is the a function; there are no outputs.

subroutine SOLDIF (ausr)

external ausr

4.2.4 SOLINP handles most of the details of input for a problem. All the arguments are outputs. Their meanings were discussed in Section 3.4.7, under USRINl.

subroutine SOLINP (oldsol, nt. nv. nc . nb, npde)

logical oldsol

integer nt . nv . nc . nb

4.2.5 SOLITR handles the nonlinear iterations; it is the workhorse routine. The inputs

are the names of the three subroutines describing the PDEs. SOLITR asks for values of

iteration parameters and then proceeds.

subroutine SOLITRCausr, fusr, gusr)

external ausr. fusr, gusr

ausr, fusr, and gusr are the names of subroutines that describe the PDEs and the

boundary conditions.

4.2.6 SOLREF is used either to start a new solution or to refine the mesh for an existing

solution.

subroutine SOLREF (ausr, fusr. gusr. nevlvl)

external ausr. fusr, gusr

ausr, fusr, and gusr are the names of subroutines that describe the PDEs and the

boundziry conditions.

newlvl should be 0 to get started on the same level, or 1 to refine to the next higher

level.

4.2.7 SOLSAV saves the solution values eind the data structure by writing them to a file in a format suitable for later use by B2DE. SOLSAV calls USRSAV in case the user wants to

write anything for later retrieval by USRGET.

subroutine SOLSAV

4.3 PLTxx Routines

27
PLTxxx routines are callable by the user, but should not be changed. In the default driver,

only the corresponding DRVxxx routines call the PLTxxx routines. The five routines are PLTCON, PLTFLW, PLTPRO, PLTSUR, and PLTTRI. They are described separtately.

4.3.1 PLTCON does a contour plot of a PDE variable. The contour values are printed to

the standard output, but do not appear on the plot. The plot uses the current color (or

color palette for filled contours) , window, and magnification.

subroutine PLTCON (nth, ntype, fill, istcol, nlev, ulev)

integer nth, ntype, istcol, nlev

real ulev(*)

logical fill

nth: Plot contours of the n-th PDE.

ntype: If ntype is 1, plot contours of u„(x, y); otherwise plot contours of | Vun(i, y) |. fill: If fill is .TRUE. , the spaces between contours are filled in; otherwise only the

contour lines lire drawn.

istcol is the starting color for filling spaces between contours, istcol is ignored if fill is .FALSE.

nlev: If nlev is positive, plot nlev contours, using automatic scaling, spacing the

contours evenly between, but not including, the maximum and minimum values of

the plot variable; ulev is ignored by PLTCON. If nlev is negative, plot -nlev contours,

using the contour values in the ulev array. If nlev is zero, just find and print the

minimum and maximum values of the plot variable.

4.3.2 PLTFLW draws flow lines for a pde variable. The amount of positive and negative

flow is printed to the standard output. The plot uses the current colors, window, and

magnification.

subroutine PLTFLW(axy, nth, npf , nflow, xyval, rp)

real xyval ()

external axy

axy is the name of the a subroutine,

nth is as described under PLTCON.

npf is the number of points defining the crossing line,

xyval is an airray of the points, first the npf z-values, then

nflow is the number of flow lines to be plotted.

4.3.3 PLTPRO draws profile plots, the outline of slices through a

uses the current color emd window, but uses magnification of 1.

subroutine PLTPRO (nth, nprof, dir)

integer nth, nprof

real dir (2)

nth is as described under PLTCON.

the npf y-values.

projection plot. The plot

28
nprof is the number of profiles.

dir are the indices of the slicing plane, as discussed in Section 3.3.4.

4.3.4 PLTSUR draws a surface projection plot of a PDE variable. The plot uses the current

color and window, but uses magnification of 1. subroutine PLTSURCnsurf , level, nth, ntype, vievrpt, ngrid,

* nlev, ulev)

integer nsurf, level, nth, ntype, ngrid (2) , nlev

real viewpt(3), ulev(*)

nsurf controls what is plotted on the surface. If nsurf is 1, trizingles of the specified

level are plotted. If nsiirf is 2, contours are plotted. If nsurf is 3, x-y grids are

plotted.

level is the level of solution to be plotted,

nth and ntype a^e as described under PLTCON.

viewpt contains the (r, y, 2) components of the view point. The surface projection

will appear as if it were at the origin, seen from this point.

nlev amd ulev are used if only nsurf is 2. Their meauings are as given under PLTCON.

4.3.5 PLTTRI draws a plot of one triangulation. The plot uses the current color, window,

and magnification.

subroutine PLTTRI (level)

integer level

level is the level of the mesh to be plotted.

4.4 Miscellaneous routines

4.4.1 lOGET is used to get the vadue of one of the input/output pau’ameters. See the USRMAI

section for a list of the meatnings.

integer function ioget (ionum)

integer ionum

ionum is the index of the paxauneter whose value is desired.

4.4.2 lOSET is used to set am input/output pairameter. See the USRMAI section for a list of the meanings.

subroutine ioset (ionum, ioval)

integer ionum. ioval

ionum is the index of the pairauneter to be set. ioval is the value to set. 4.4.3 FREEIN is a subroutine to do free-format input. The calling prograim asks for some

real numbers amd/or integers to be read from am input file. FREEIN reads a line at a time,

parses it, evaluates the numbers, and does error checking. If an error occurs and the session

is interactive, FREEIN prompts the user for a new line. Items are sepau-ated by spaces; lines aue read until all the requested items are obtained. If there are more items on a line than

29
the call to FREEIN requires, the excess is ignored. A * or # character majks the beginning

of a comment; the character and the rest of the line are ignored.

subroutine FREEIN (ioread, ioerr, nitems, itype, ibuff , rbuff)

integer ioread, ioerr, nitems, itype (nitems) , ibuff (nitems)

real rbuff (nitems)

ioread is the file number from which to read the data,

ioerr is the file number on which to write error messages,

nitems is the number of items to be read.

itype(j) is the type of the j-th item, 1 for an integer and 2 for a real number.

ibuff and rbuff are arrays into which the data is put. If itype (j ) is 1, the j-th item

is interpreted as an integer and put into ibuff (j); if itype(j) is 2, the j-th item is interpreted as a real number and put into rbuff ( j )

4.4.4 FLTGET is provided for convenience; it returns one real value typed at the terminal,

using FREEIN.

real function fltget(junk)

j unk is a dummy parameter.

4.4.5 INTGET is provided for convenience; it returns one integer value typed at the terminal,

using FREEIN.

integer function intget(junk)

junk is a dummy parameter.

4.4.6 EVAL evaluates PDE variables for passing to other programs.

subroutine EVAL (kind. nth. npde, npts, x. y. u. dudx, dudy)

real x(npts) , y(npts) , u(npts) , dudx(npts) , dudy(npts)

If kind is 1, only the PDE variable is evaluated and returned in u; dudx zmd dudy are not referenced. If kind is 2, only the gradient of the PDE variable is evaluated and

retiirned in dudx and dudy; u is not referenced. If kind is 3, both the PDE variable

and its gradient are evaluated.

The nth PDE variable is evaluated.

npde is the number of PDE variables.

npts is the number of points at which to evaluted the nth PDE.

X and y are the points at which to evaluate the nth PDE.

If kind is 1 or 3, the PDE values are returned in u.

If kind is 2 or 3, the PDE gradient values zu'e returned in dudx and dudy.

5. Installing B2DE on Your Computer

Most of B2DE is standard Fortran 77; it should not need changes by the user euid should

not be changed. All known machine dependencies have been isolated for easy changing.

30
5.1

The PORT Stack

Working space is managed inside B2DE using the PORT stack [Fox78]. The size of the

stack determines how large a problem can be solved. The stack takes advantage of a

feature which is in most Fortrans, but which is not in the Fortran 77 standard: the size of a labelled common block is determined by the size in the main program, or by the size in

the first program loaded. To change the size of the stack, only the main program needs to

be edited and re-compiled. The following block of code occurs many times within B2DE,

defining the stack common, cstak, to be a nominal size.

parameter (irsize=1000)

parameter (idsize- 500)

common/cstak/dstak

double precision dstak(idsize)

real r(irsize)

integer i(irsize)

logical IstedcCirsize)

equivalence (dstak(l) .r(l)) . (dstak(l) .i(D) , (dstak(l) .Istak(l)

)

save /cstedc/

Your main program will define the true size of the stack.

5.2 Machine-Dependent Programs

To be as portable as possible, B2DE uses the machine constants from the PORT library

[Fox78j. The machine constant functions IlMACH and RIMACH must be changed as appropriate for your computer. The distributed version is correct for the DEC VAX 11/700 series

and the VAX/VMS operating system. Constants for many other computers are in the

routines, but commented out. See [Fox78] for instructions on changing machine constants

for computers not listed.

B2DE \ises two timing routines. Each should calculate an elapsed amoimt of time used in

milliseconds. ISTIME is called at the beginning of an activity to be timed, and IFTIME is called at the conclusion. The difference is used as the elapsed time.

integer function iftime ( j\ink ) junk is a dummy parameter.

integer f\mction istime ( junk ) junk is a dummy parameter.

5.3 Terminal-Dependent Programs

TRMxxx routines are terminal-dependent. The default versions work for DEC VT125 ter- minals iind Hewlett Packard 264x and 262x terminals. New versions must be written for

other terminals. Some TRMxxx routines are not called if plotting to a file is being done, and

therefore are not necessary if only batch operation is planned, although dummy versions

must be provided to satisfy most operating systems.

5.2.1 TRMCOL should set the color or the line type for the plots. If the terminal does not

support color or alternate line types, TRMCOL need not do anything. TRMCOL is called by

SETCOL, which is called by DRVCOL; DRVCOL only allows non-negative values for the color.

31
subroutine TRMCOL(newcol)

integer newcol

newcol is the new color or line type.

5.2.2 TRMDRW should draw a line from the last point reached by TRMDRW or TRMMOV to the

point (ix, iy). The point (ix, iy) is guaranteed to be within the limits set by TRMSCL.

TRMDRW is not called if plotting to a file is being done. The line should be drawn in the

current color or line type as last set by TRMCOL.

subroutine TRMDRW (ix, iy)

integer ix, iy

ix and iy are the screen coordinates of the new point.

5.2.3 TRMMOV should move the current point to the point (ix, iy). No line should be

drawn. The point (ix, iy) is guaranteed to be within the limits set by TRMSCL. TRMMOV

is not called if plotting to a file is being done.

subroutine TRMMOV (ix, iy)

integer ix, iy

ix and iy eire the screen coordinates of the new point.

5.2.4 TRMPAL changes the color palette available to the user. It has no arguments,

subroutine TRMPAL

5.2.5 TRMPLT is a plotting utility. Its function depends upon its argument, as described

below.

subroutine TRMPLT (if lag)

integer if lag

if lag = 1: TRMPLT should do whatever is necessary to initialize all plotting for the

session.

if lag = 2: TRMPLT should do whatever is necessary to begin a plot, amd should erase the screen.

if lag = 3: TRMPLT should do whatever is necessary to begin a plot, but should not

erase the screen. This is not possible on some terminals.

if lag = 4: TRMPLT should do whatever is necessary to finish a plot. If the session is interactive, and if plots are done on the same screen as text, TRMPLT should wait for a signal from the user before returning.

iflag = 5: TRMPLT should do whatever is necessary to finish aill plotting for the

session.

5.2.6 TRMSCL should give the minimum and maximum screen coordinates for terminal type

iterm. If iterm is zero or am unknown type, the values for plotting to a file should be

returned.

integer function TRMSCL (iterm. iflag)

integer iterm. iflag

iflag = 1: TRMSCL should return the minimum x value.

32
if lag = 2: TRMSCL should return the maximum x value,

if lag = 3: TRMSCL should return the minimum y value,

if lag = 4: TRMSCL should return the maximum y value.

6. Examples

6.1 Laplace’s Equation

The equation is Laplace’s equation, the simplest elliptic PDE. The region is a square with

two comers roimded off, as in Fig. 1.

V • Vu = 0

The boimdary conditions 2ire u = 1 on edges 1 zmd 2, u = 0 on edge 5 , zmd dujdv — Q on

sides 3 , 4, and 6. The solution is singular at vertices 1 and 6.

Appendix A is a program to solve this example; Appendix B is an input data file suitable

for use with the default driver; and Appendix C is a sample run. A typical initiad triangulation supplied by a user is shown in Figure 2. Sides 1 and 4 are defined as circular arcs, but the curvature does not show up imtil the sides are refined.

Figure 3a is the actual coarsest triangulation used; the input data file specifies initial

refinement arotmd the singularities. Figures 3b amd 3c are refined meshes generated by

B2DE in the sample run. Note that the adaptive refinement concentrated new triangles neax the singularities. Figure 4 shows other plots froni the sample run: 4a is a contour

plot; 4b is a fiow line plot; and 4c is a profile plot.

6.2 A Nonlinear Example

This is a concocted set of three nonlinear elliptic PDEs. The solutions are smooth, and the

adaptive mesh generates uniform triangulations. The region is the rectangle — 1 < a: < 1 , 0 < y < 1.

V • (VU2) = 0

V • (U1U3VU3) = 0

The boundary conditions are all Dirichlet:

U1U3U3 = 1

U3 = e*

U3 =

33
The solution is

U2 = e®

Us =

Appendix D is a program to solve this example; Appendix E is an input data file suitable

for use with the default driver; and Appendix F is a sample run. A typical initial triangulation supplied by a user is shown in Figure 5a. Figure 5b is the

actual coarsest triangulation used; the input data file specifies initial refinement around the

left comer. (There is no singularity in the problem; this initial refinement is unnecessary.)

Figures 5c zuid 5d are refined meshes generated by B2DE in the sample run. Note that the

adaptive refinement made the mesh uniform; the incorrect initial mesh has been overcome.

Figure 6 shows other plots from the sample nm: 6a is a surface plot showing triangles; 6b

is a surfax:e plot showing contours, with hidden lines not removed; and 6c is a surface plot

showing grid lines.

References

[Bank79] R. E. Bank and A. H. Sherman, “PLTMG Users’ Guide,” CNA 152, Center for

Numerical Analysis, The University of Texas at Austin, September, 1979.

[Bank82| R. E. Bank, “PLTMG Users’ Guide, Jime, 1981 version,” Technical Report,

Department of Mathematics, University of California at San Diego, Augiist, 1982.

[BlueSS] J. L. Blue and C. L. Wilson, “Two-Dimensional Analysis of Semiconductor Devices

Using General-Purpose Interactive PDE Software,” IEEE Trans. Electron Devices, vol.

ED-30, pp. 1056-1070, 1983.

[Bran77] A. Brandt, “Multi-Level Adaptive Solutions to Boundary-Value Problems,” Math.

Comput., vol. 31, pp. 333-390, 1977.

[Fox78| P. A. Fox, A. D. Hall, md N. L. Schryer, “The PORT Mathematical Subroutine

Library,” ACM Trans. Mathematical Software, vol. 4, pp. 104-126, 1978.

[Kell86| E. F. Kelley, R. E. Hebner, W. E. Anderson, J. A. Lechner, and J. L. Blue, “The

Effect of an Oil-Paper Interface Parallel to an Electric Field on the Breakdown Voltage at

Elevated Temperatures.” Submitted to IEEE Trans. Power Apparatus and Systems.

[Kern75] B. W. Kemighan, “RATFOR- A Preprocessor for a Rational Fortran,” Software - Practice and Experience, vol. 5, pp. 395-406, 1975.

[Sher78| A. H. Sherman, “Algorithms for Sparse Gaussim Elimination with Partiad Pivoting,” ACM Trans. Mathematical Software, vol. 4, pp. 330- 338, 1978.

[Stra73] G. Strang emd G. J. Fix, An Analysis of the Finite Element Method (Prentice-Hall,

Englewood Cliffs, New Jersey, 1973).

[Wils85] C. L. Wilson, P. Roitman, and J. L. Blue, “High Accuracy Physical Modeling of

Submicron MOSFETs,” IEEE Trans. Electron Devices, vol. ED-32, pp. 1246-1258, 1985.

Figure Captions

34
1. Region for Laplace’s equation example. The boundary conditions are u = 1 on edges 1 and 2, u = 0 on edge 5, and duifdv = 0 on sides 3, 4, and 6.

2. User-supplied triangulation of region. The curvature of sides 1 and 4 will show up in

refinements of the region.

3a. Initial triangulation of region. Vertices 1 and 6 were specified as refinement level 3; other vertices were level 1. 3b. Triangulation after first adaptive refinement.

3c. Triangulation after second adaptive refinement.

4a. Contours of the approximate solution on the mesh of Fig. 3c. The contours ajre imiformly

spaced.

4b. Flow line plot of the approximate solution on the mesh of Fig. 3c. Equal fiux fiows

between each pair of fiow lines.

4c. Profile plot of the approximate solution on the mesh of Fig. 3c, with slicing direction

1 0. 5a. User-supplied triangulation of region for example 2. 5b. Initial triangulation of region. Vertex 4 was specified as refinement level 3; other vertices

were level 1. 5c. Triangulation edter first adaptive refinement.

5d. Triangulation after second adaptive refinement.

6a. Surface plot of error in PDE 2, on mesh of Fig. 5d. Triangles are shown, ajid hidden

lines are removed.

6b. Surface plot of error in PDE 2, on mesh of Fig. 5d. Contours are shown, and hidden

lines are not removed.

6c. Surfax:e plot of error in PDE 2, on mesh of Fig. 5d. Grid lines are shown, 21 along x and 11 along y, and hidden lines are removed.

35
Appendix A. Program to Solve Laplace’s Equation Example

c

c main program to drive nonlinear multi-variable code to solve

c Laplace’s equation, adapted from umain.for

c

c the user-supplied routines that determine the PDEs

c external alap, flap, glap

c

c i/o value array

c integer io(9)

cc the stack must be declared and initialized in the main program,

c it is ’save’ed for safety’s sed^e. the equivalence is required on

c some con^uters to get the word alignment correct,

c parameter (irsize = 100000)

parameter (idsize = 60000)

common/cstedc/dstedc

double precision dstak(idsize)

real rstak(irsize)

equivalence (dstedcCl) . rstak(l))

save /cstak/

cc initialize the stack to tell it the true total size

c call istkinCirsize, 3)

c

c specify i/o values

c

c io(l) = file number for reading from the terminal

c io(2) » file number for writing to the terminal

c io(3) * terminal type, used by graphics programs

c io(l) = 5

io(2) * 6

io(3) » 126

cc io(4) = file number for reading input data;

c only used by usrinl . usrin2 , and usrin3

c io(6) = file number for writing plot data, if not to terminal

c io(6) * file number for writing print data, if not to terminal

c io(4) = 20

io(6) * 21

io(6) = 22

c

36
c io(7) = file number for reading old solution data

c io(8) = file number for writing new solution data

c io(9) = file number for writing scratch data

c For reading aind writing solution data, the data is written as

c unformatted, or binary, if the io value is positive, and as

c formatted if the io value is negative,

c io(7) = -23

io(8) = -24

io(9) = 25

cc install the i/o values

c do 10 j = 1. 9

call iosetCj, io(j))

10 continue

c

c open data files,

c openCunit = io(4), file = ’inputl.dat’, status = ’old’)

openCunit = io(5) , file = 'plot.daf , status = ’unknown’)

openCunit * io(6), file = ’print.dat’ , status = ’unknown’)

c if ( io(7) .gt. 0) then

openCunit = ioC7), file » ’insol.dat’, status = ’unknown’,

* form* ’unformatted’)

else

openCunit = -ioC7) , file = ’insol.dat’, status = ’unknown’,

* form* ’formatted’)

endif

c if C ioC8) .gt. 0) then

openCunit * ioC8) , file * ’outsol.dat’, status = ’unknown’,

* form* ’unformatted’)

else

openCunit = -ioC8), file = ’outsol.dat’, status = ’unknown’,

* form* ’formatted*)

endif

c openCunit * ioC9) , file * ’scrtch.dat’, status * ’\mknown’,

* form* ’unformatted’)

c

c teike it away . . . c call drvallCalap, flap, glap)

c

c close the files, just to be safe

c

37
close (io(4))

close(io(5))

close(io(6))

close (labs (io (7)))

close(iabs(io(8) ))

close (io (9))

c stop

end

c

c subroutine alap(npde, ndim, dopde, x, y, i, u, dudx, dudy,

* a, dadu, dadux, daduy)

real x, y, u(*) , dudx(*) , dudy(*)

* a(*)

daduCndim,*) , daduxCndim,*) , daduy(ndim, )

logical dopde (npde)

c

c inputs:

c npde xthe number of pdes in the problem

c ndim the first dimension of arrays dadu. dadux, eind daduy

c dopde dopde (n) is .TRUE, if the nth is active,

c 7 the coordinates of the point at which to evaluate

c i the macro triangle number that the point is in

c u the current pde values at the point

c dudx. dudy the gradient of the pdes at the point

c

c outputs

c a the values of the a coefficient

c dadu the first derivatives of a with respect to u

c dadux the first derivatives of a with respect to du/dx

c daduy the first derivatives of a with respect to du/dy

c a(l) * 1.

daduCl, 1) * O.eO

daduxCl, 1) » O.eO

daduy (1, 1) =* O.eO

c return

end

cc subroutine flap (npde, ndim, dopde, x, y, i, u, dudx, dudy,

f, dfdu, dfd\ix, dfduy)

real x, y, u(*)

dudx(*) , dudy(*)

f(*). dfdu(ndim, *) , df dux (ndim, *) , dfduy (ndim, •)

logical dopde (npde)

cc

38
c inputs:

c npde xthe number of pdes in the problem

c ndim the first dimension of airrays dadu, dadux, and daduy

c dopde dopde(n) is .TRUE, if the nth is active,

c 3C, y the coordinates of the point at which to evaluate

c i the macro triaoigle number that the point is in

c u the current pde values at the point

c dudx.dudy the gradient of the pdes at the point

c

c outputs:

c f the values of the f coefficient

c dfdu the first derivatives of f with respect to u

c dfdux the first derivatives of f with respect to du/dx

c dfduy the first derivatives of f with respect to du/dy

c f(l) = O.eO

dfdu (1. 1) = O.eO

dfduxCl, 1) = O.eO

dfduyCl, 1) » O.eO

c return

end

c subroutine glapCnpde, ndim, dopde. x. y, i, j, u. g, dgdu)

real x, y, u(*) , g(*) . dgduCndim,*)

logical dopde (npde)

c

c inputs:

c npde

c ndim

c dopde

c X. y

c i

c j c u

c

c outputs

c g . c dgdu

cc

c square with 2 round comers

c if (i.eq.l .or. i.eq.2) then

g(l) = u(l) - 1.0

dgdud, 1) » 1.0

else if (i .eq. 5) then

g(l) » u(l)

xthe number of pdes in the problem

the first dimension of arrays dadu, dadux, and daduy

dopde (n) is .TRUE, if the nth is active,

the coordinates of the point at which to evaluate

the macro triangle number that the point is in

the side of the macro triangle that the point is on

the current pde values at the point

the values of the g coefficient

the first derivatives of g with respect to u

39
dgduCl. 1) = 1.0

else

g(l) = 0.0

dgduCl, 1) =0.0

endif

return

end

40
Oio

Appendix B. File “inputl.dat” for Laplace’s Equation Example

# iziput.dat for Laplace example

It square with two rounded comers

# old solution flag

7 2 6 1 # nt nv nc nb npde

# vert Ivl X y

1 3 -0.5 0.0

2 1 0.0 0.5

3 1 0.5 0.5

4 1 0.0 0.0

5 1 0.5 0.0

6 3 0.0 -0.5

7 1 -0.5 -0.5

# curved edges

1 .35355339 -.35355339

2 -.35355339 .35355339

# triangles

112 4

2 2 3 4

3 3 4 5

4 4 5 6

5 4 6 7

6 14 7 . 11321 # be tri edge curv be

2 2 3 0 1 3 3 2 0 2

4 4 112

5 5 10 1 6 6 2 0 2

2 8 0 # maxm maxr ittype

4 17 It Ivlref Ivlusr maxi

-30 # ithresh

2 # iprsw

80000 # work space

41
Appendix C. Transcript of interactive session for Laplace’s Equation Example

Output from the computer is in typewriter font. Numbers typed by the user are in

boldface and have a • at the beginning of the line.

refsol entered at level 1 7 vertices

grid 0.040 sec. 26 vertices

asmbnl 0.030 sec. l=print, 2=plot, 3=iter, 4=refine, 5=comp, 6=save, 7=diff, 8=end, 9=change

• 3

enter niter, ieps. relerr, abserr, jsmax, paraun

• 1 0 0 0 1 0

Iter 1

1 residual L2 norm 4.397E-01

factor 0.040 sec. 26 eqs.

asmbnl 0 . 250 sec

mg 0.000 sec. GE

1 newton norm, solution norm, relerr l.OOOE+00 O.OOOE+00 l.OOOE+00

0=continue, l=print, 2=plot

• 0

0 t 0.000000 L2(rhs) 4.397E-01 0.110 sec.

1 t 1.000000 L2(rhs) 1.464E-07 0.110 sec. l»print, 2=plot, 3=iter, 4=refine, 5=comp, 6=save, 7=diff, 8=end, 9=change

• 4

refsol re-entered at level 2 26 vertices

asmbnl 0 . 030 sec

1 estimated maximum error 8.01E-03

estimated individual PDE errors

0=continue , l=print , 2*plot

• 0

level 1 pde 1 \morm, enorm 7.06E-01

min error / max error 1.44 percent

4 of 30 triangles ( 13.33) above 50.00 percent

9 of 30 triangles ( 30.00)

8 of 30 triamgles ( 26.67)

refine 0.300 sec. level 2

l=print, 2=plot, 3=iter, 4=refine,

• 5

2.86E-02 eat. digits 1.39

above 25.00 percent

above 37.60 percent

52 vertices

5=comp, 6=save , 7=diff, 8»end, 9=change

compressed to level 1 l=print, 2=plot, 3=iter, 4=refine, 5=comp, 6=save, 7=dif f , 8=end, 9»change

• 3

enter niter, ieps, relerr, abserr, jsmax. param

• 100010

Iter 1

1 residual L2 norm 3.650E-02

factor 0.080 sec. 52 eqs.

asmbnl 0 . 480 sec

mg 0.000 sec. GE

1 newton norm. solution norm, relerr 7.664E-02 l.OOOE^OO 7.564E-02

0=continue, l=print, 2=plot

42
• 0

0 t 0.000000 L2(rhs) 3.650E-02 0.260 sec.

1 t 1.000000 L2(rhs) 7.992E-08 0.240 sec. l=print, 2=plot, 3=iter, 4=refine, 6=comp, 6=save, 7=diff, 8=end, 9=ch2mge

• 4

ref sol re-entered at level 2 52 vertices

asmbnl 0.050 sec.

1 estimated maximum error 2.64E-03

estimated individual PDE errors

0=continue, l=print, 2=plot

• 0

level 1 pde 1 min error / max

9 of 72

26 of 72

14 of 72

19 of 72

23 of 72

20 of 72

21 of 72

22 of 72

refine 0.720

l=print , 2=plot

• 3

unorm, enorm 7.04E-01

error 5.91 percent

triangles ( 12.50) above 50.00 percent

1.04E-02 est. digits 1.83

triangles ( 36.11)

triangles ( 19.44)

triangles ( 26.39)

triangles ( 31.94)

triangles ( 27.78)

trieoigles ( 29.17)

triangles ( 30.56)

sec . level 2

, 3=iter, 4=refine,

above 25.00 percent

above 37.50 percent

above 31.25 percent

above 28.13 percent

above 29.69 percent

above 28.91 percent

above 28.52 percent

103 vertices

5=comp, 6=save, 7=diff, 8=end, 9=change

enter niter, ieps, relerr, abserr, jsmax, param

• 1 10 0 0 1 0

Iter 1 1 residual L2 norm 1.761E-02

factor 0.080 sec . 52 eqs

asmbnl 1.240 sec

mg 2 7.619E-01 3.472E-02 4.557E-02

mg 2 3.472E-02 4.404E-03 1.268E-01

mg 2 4.404E-03 5.974E-04 1.357E-01

mg 0.160 sec. GS reduction 7.84E-04

1 newton norm, solution norm, relerr 5.398E-02 l.OOOE+00 5.398E-02

0=*continue, l=print, 2*plot

• 0

0 t 0.000000 L2(rhs) 1.761E-02 0.510 sec.

1 t 1.000000 L2(rhs) 1.659E-06 0.500 sec. l=print, 2*plot, 3=iter, 4*refine, 5=comp, 6=save, 7=diff, 8=end, 9=change

• 1

l=8ummary , 2»timing, 3=storage, 4=vertices, 5=pde variables, 6=settings

• 1

n Ivl vertices triangles unorm enorm digits 11 52 114 7.04E-01 1.04E-02 1.83

1 2 103 258

l=summary, 2=timing, 3=storage, 4=vertices, 5=pde variables, 6=settings

• 2

43
timing distribution

plot

seconds

0.000

assembly 2.080

factor 0.200

mg 0.160

rhs check 1.730

refine 1.060

total 5.230

l=summary2=timing

• 3

103 vertices

2094 vertices

263 trieuigles

6583 trieingles

l=3ummary2=timing

• 6

percent

0.00

39.77

3.82

3.06

33.08

20.27

100.00

3=storage, 4=vertices,

used.

allowed.

used.

allowed.

3=storage, 4=vertices,

5=pde variables , 6=setting3

6=pde variables, 6=setting3

l=color, 2=erase, 3=magnify, 4=window, 5=file print/plot, 6=transform

• 3

enter magnification

• 3

enter x-center and y-center (0 to 100)

• 50 0

l=color, 2=erase, 3=magnify, 4=window, 6=file print/plot, 6=transform

• .0 l«a\nranary, 2=timing, 3=storage, 4=vertices, 6=pde variables, 6=“settings

• 5

l=on grid. 2=at vertices, 3=on grid as scaled integers

• 3

enter pde no., l=value only, 2=gradient only, 3= both

• 11

nx, ny

• 77

X limits -1.66667E-01 1.66667E-01

y limits -6.00000E-01 -1.66667E-01

pde 1 scaled values from -6.086E-08 to 6.492E-01 scale 5.492E-01

90 90 91 93 96 97 100

78 78 80 82 86 89 92

64 66 67 70 76 80 86

49 60 54 68 66 71 77

34 36 39 44 64 62 70

17 19 21 28 42 66 66

0 0 0 0 9999 9999 9999

l=on

• 0

grid. 2=at vertices

3“on grid as scaled integers

l=summary

• 0

2-timing, 3=storage , 4=vertices, 6=pde variables.

44
l=print, 2=plot, 3=iter, 4=refine, 5=comp, 6=save, 7=diff, 8=end, 9=change

• 2

I=tri2uigle8 , 2=conto\irs, 3=surface, 4=profiles, 6=flow lines, 6=settings

• 6

l=color, 2=erase, 3=magnify, 4=window, 6=file print/plot, 6=transform

• 3

enter magnification

• 1

l*color, 2=erase, 3=magnify, 4=window, 6=file print /plot, 6=transform

• 0

l=triangles, 2=contours, 3=surface, 4=profiles, 5=flow lines, 6=settings

• 1

enter level

• 1

plot 0 . 360 sec

enter level

• 2

plot 0.620 sec. enter level

• 0

l=triangles, 2=contours, 3=surface, 4=profiles, 5=flow lines, 6=settings

• 2

pde no., l=function/2='abs (grad) , l=fill/2=no fill

112

number of contours

• 11

min and max values O.OOOOOOE+00 l.OOOOOOE+00

contour values

8.333E-02 1.667E-01 2.600E-01 3.333E-01 4.167E-01

6.000E-01 6.833E-01 6.667E-01 7.600E-01 8.333E-01

9.167E-01

plot 0 . 480 sec

pde no., l=function/2=abs(grad) , l=fill/2=no fill

• 0 0 0

l=triangles, 2=contours, 3=surface, 4=profiles, 5=flow lines, 6=settings

• 3

pde no. , l=function/2=abs(grad)

• 00

l=triangles, 2=contours, 3=8tirfac8, 4=profiles, 6=flow lines, 6=settings

• 5

pde no., number of points, number of flow lines

• 1 2 11

enter 2 x-values followed by 2 y-values

• "•S "•S ,6

total +, - flow 9.743972E-01 O.OOOOOOE+00

plot 0.760 sec. pde no., number of points, number of flow lines

• 0 0 0

45
l=triangles , 2=contours, 3=surface, 4=profiles, 6=flow lines, 6=setting3

• 4

pde no., number of profiles

• 1 21

indices of profile planes

• 10

plot 0.910 sec. pde no., number of profiles

• 00

l=triangles, 2=contours, 3=surface, 4=profiles, 5=flow lines, 6=settings

• 0

l=print, 2=plot, 3=iter, 4=refine, 5=comp, 6=save, 7=diff, 8=end, 9=change

• 8

46
Appendix D. Program to Solve Nonlinear System Example

external aex, fex, gex

integer io(9)

parameter (irsize = 400000)

parameter (idsize = 200000)

c ommon/c Stak/dsteik

double precision dstak(idsize)

real rstak(irsize)

equivalence (dstak(l) , rstak(l))

save /cstak/

call istkinCirsize . 3)

io(l) » 5

io(2) = 6

io(3) » 125

io(4) * 20

io(5) * 21

io(6) » 22

io(7) » 23

io(8) = 24

io(9) = 26

do 10 j * 1, 9

call iosetCj, io(j))

10 continue

openCunit = io(4) , file = ’input2.dat’, status = ’old*)

openCunit * io(5), file = ’plot.dat’ , status = ’unknown’)

openCunit = io(6) , file = ’print.dat’ , status = ’unknown’)

if ( io(7) .gt. 0) then

openCunit = ioC7) , file = ’insol.dat’, status = ’unknown’,

form* ’unformatted’)

else

openCunit * -ioC7) , file = ’insol.dat’, status = ’unknown’,

form* ’formatted’)

endif

if C ioC8) .gt. 0) then

openCunit * ioC8) , file = ’outsol.dat’, status = ’unknown’,

form* ’unformatted’)

else

openCunit * -ioC8) , file = ’outsol.dat’, status = ’unknown’,

form* ’formatted’)

endif

openCunit * ioC9) , file * ’scrtch.dat’, status *

’ unknown*

form* ’unformatted’)

call drvallCaex. fex. gex)

closeCioC4)

)

closeCioCS))

closeCioCd)

close CiabsCioC7)))

closeCiabsCioC8) )

)

47
close (io (9)

)

stop

end

subroutine aexCnpde, ndim, dopde, x, y. i, u, dudx, dudy,

* a. dadu, dadux. daduy)

real x, y, u(*) , dudx(*) , dudy(*)

a(*) , daduCndim,*) , daduxCndim,*) , daduy (ndim, *)

logical dopde (npde)

do 10 j

=» 1 , npde

do 10 k = 1 , npde

dadu (j, k) = 0.0

daduxCj , k) =0.0

daduy ( j , k) = 0.0

10 continue

a(l) = 1.0 + dudx(2)**2 + dudy(2)**2 + dudx(3)**2 + dudy(3)**2

daduxCl, 2) = 2.0*dudx(2)

daduy (1, 2) = 2.0*dudy(2)

daduxCl, 3) = 2.0*dudx(3)

daduy (1. 3) = 2.0*dudy(3)

a(2) = 1.0

a(3) » u(l)*u(2)

dadu(3, 1) = u(2)

dadu(3, 2) = u(l)

ret\im

end

subroutine f ex (npde, ndim, dopde, x, y, i, u, dudx, dudy.

f, dfdu, dfdiox, dfduy)

real x, y, u(*) , dudx(*) , dudy(*)

* f(*), dfdu(ndim, *) , df dux (ndim, *) , dfduy (ndim, *)

logical dopde (npde)

do 10 j =1, npde

do 10 k = 1 . npde

dfdu (j, k) = 0.0

dfdux(j , k) =0.0

dfduy (j , k) =0.0

10 continue

f(l) = 2.0*u(l)

dfdud, 1) * 2.0

f(2) = u(2)**2 * exp(-x)

dfdu(2, 2) = 2.0 u(2) exp(-x)

f(3) = 0.0

return

end

subroutine gex(npde, ndim, dopde, x, y, i, j, u, g, dgdu)

real x, y, u(*) , g(*)

dgdu(ndim,*)

logical dopde (npde)

g(l) = u(l)*u(2)*u(3) - 1.0

dgdud, 1) = u(2)*u(3)

48
dgdu(l, 2) = u(l)*u(3)

dgdu(l, 3) = u(l)>«'u(2)

g(2) * u(2) - exp(x)

dgdu(2, 2) = 1.0

g(3) = u(3) - exp(y)

dgdu(3, 3) = 1.0

return

end

subroutine usrgssCnpde, x. y, i, u)

real u(*)

do 10 n = 1, npde

u(n) =1.0

10 continue

return

end

subroutine usrshf (ixtype , npde, nv, u, xmew, vx, vy)

real u(npde,nv), unew(nv) , vx(nv) , vy(nv)

if (ixtype .le. npde) then

nth = ixtype

do 10 i = 1, nv

unew(i) = u(nth, i) - utrue(nth, vx(i) , vy(i))

10 continue

else

do 20 i = 1 , nv

unew(i) = u(l,i) u(2,i) * u(3,i)

20 continue

endif

return

end

real function utrueCnth, x, y)

real x, y

if (nth .eq. 1) utrue = exp(-x-y)

if (nth .eq. 2) utrue » exp(x)

if (nth .eq. 3) utrue » exp(y)

return

end

49
Appendix E. File “input2.dat’* for Nonlinear System Example

It two squares

It dirichlet b.c. on all sides

0 It old solution flag

It nt nv nc nb npde

4 6 0 6 3

# vertices

# i Ivl X y

1 1 - 1.0 1.0

2 1 0.0 1.0

3 1 1.0 1.0

4 3 -1.0 0.0

6 1 0.0 0.0

6 1 1.0 0.0

# curved edge centers

# (none)

# triangles

# i vl v2 v3

112 4

2 2 3 5

3 2 4 5

4 3 5 6

# boundary conditions:

# i triangle edge curve b.c.

112 0 111

2 13 0 111

3 2 3 0 1 1 1 4 4 2 0 1 1 1 5 4 10 111

6 3 10 111

2 8 0 # maxm maxr ittype

407 # Ivlref Ivlusr Ivlmax

20 # threshold percentage

1 # print level

350000 # work space

50
Appendix F. Transcript of interactive session for Nonlinear System Example

Output from the computer is in typewriter font. Numbers typed by the user are in

boldface and have a • at the beginning of the line.

refsol entered at level 1 6 vertices

grid 16 vertices

l=print, 2=plot, 3=iter, 4=refine, 5=comp, 6=save , 7=diff, 8=end, 9=change

• 3

enter niter, ieps, relerr, abserr, jsmaoc, param

• 10 0 0 .1 5 0

Iter 1

1 newton nom. solution norm, relerr 3 . 437E+00 l.OOOE+00 3.437E+00

2 newton norm. solution norm, relerr 1.718E+00 l.OOOE+00 1.718E+00

3 newton norm. solution norm, relerr

using newton step fraction 2.500E1.718E+00

•01

l.OOOE+00 1.718E+00

Iter 2

1 newton norm, solution norm, relerr 6.817E-01 1 . 158E+00 5.887E-01

2 newton norm, solution nom, relerr 1.289E+00 1.430E+00 9.015E-01

3 newton norm. solution nom, relerr 1.289E+00 1.430E+00 9.015E-01

Iter 3

1 newton norm, solution nom, relerr 8.786E-01 1 . 840E+00 4.775E-01

2 newton norm, solution nom, relerr 8.669E-03 2.718E+00 3.189E-03

3 newton

Iter 4

norm. solution nom. relerr 5.293E-02 2.718E+00 1.947E-02

1 newton norm, solution nom, relerr 6.897E-03 2.718E+00 2.169E-03

2 newton norm, solution nom, relerr 4.080E-06 2.718E+00 1.501E-06

3 newton norm, solution nom, relerr 1.173E-02 2.718E+00 4.315E-03

l=print, 2=plot, 3=iter, 4=refine, 5=conqp, 6=save, 7=diff, 8=end, 9=change

• 4

refsol re-entered at level

1 estimated maximum error

2 estimated maximum error

3 estimated maximum error

2 16 vertices

7.12E-02

6.13E-02

6.38E-02

estimated individual PDE errors

0=continue, l=print, 2=plot

• 0

estimated total error (as PDE #1)

0=continue , l=print , 2»plot

• 0

refine to level 2 27 vertices

l=print, 2=plot, 3=iter, 4=refine, 5=comp, 6=save, 7=diff , 8=end, 9=change

• 5

compressed to level 1 l=print, 2=plot, 3=iter, 4=refine, 6=comp, 6=save, 7=diff , 8=end, 9=change

• 3

enter niter, ieps. relerr. abserr, jsmax. param

• 10 0 0 .01 5 0

Iter 1

1 newton norm, solution norm, relerr 2.014E-01 2.718E+00 7.411E-02

51
2 newton norm, solution norm, relerr 2.14©E-©1 2.718E+00 7.874E-©2

3 newton

Iter 2

norm, solution norm, relerr 2.1©4E-©1 2.718E+©© 7.741E-©2

1 newton norm, solution norm, relerr l.©©7E-©2 2.718E+©© 3.7©6E-©3

2 newton norm, solution norm, relerr 9.477E-©4 2.718E+©© 3.486E-©4

3 newton

Iter 3

norm, solution norm, relerr 3.569E-©2 2.718E+©© 1.313E-©2

1 newton norm, solution norm, relerr 7 . 878E-©6 2.718E+©© 2.898E-©6

2 newton norm, solution norm, relerr 4.161E-©8 2.718E+©© 1.527E-©8

3 newton norm, solution norm, relerr 1.327E-©4 2.718E+©© 4.883E-©5

l=print, 2=plot, 3=iter, 4=refine, 5=comp, 6=save, 7=diff, 8=end. 9=change

• 4

ref sol re-entered at level

1 estimated maximum error

2 estimated maximiim error

3 estimated meacimum error

2 27 vertices

1.70E-02

6.18E-02

7.46E-02

estimated individual PDE errors

©“continue , 1 “print , 2=plot

• 0

estimated total error (as PDE #1)

©“continue , 1 “print , 2=plot

• 0

refine to level 2 45 vertices

Imprint, 2*plot, 3“iter, 4=refine, 6=comp, 6=save , 7»diff, 8=*end, 9*change

• 3

enter niter, ieps, relerr, abserr, jsmax, param

• 10 10 0 .001 5 0

Iter 1 Chkmtx 136 equations. Worst equation 1© 9.894E-©2

This set of equations may be unsuitable for Gauss-Seidel iteration.

1 newton norm, solution norm, relerr 7.868E-©2 2.718E+©© 2.891E-©2

2 newton norm, solution norm, relerr 6.731E-©2 2.718E+©© 2.476E-©2

3 nevrton norm, solution norm, relerr 6.66©E-©2 2.718E-*-©© 2.446E-©2

Iter 2

Chkmtx 135 equations. Worst equation 1© 9.894E-©2

This set of equations may be unsuitable for Gauss-Seidel iteration.

1 newton norm, solution norm, relerr 1.728E-©3 2.718E+©© 6.357E-©4

2 newton norm, solution norm, relerr 9.©14E-©6 2.718E-*-©© 3.316E-©6

3 newton norm, solution norm, relerr 4.439E-©3 2.718E-*-©© 1.633E-©3

Iter 3

Chkmtx 135 equations. Worst equation 1© 9.894E-©2

This set of equations may be unsuitable for Gauss-Seidel iteration.

1 newton norm, solution norm, relerr 6.482E-©6 2.718E+©© 2.385E-©6

2 newton norm, solution norm, relerr 1.491E-©7 2.718E+©© 5.485E-©8

3 newton norm, solution norm, relerr 1.462E-©6 2.718E-*-©© 5.343E-06

l“print, 2“plot, 3=iter, 4“refine, 6“comp, 6“save , 7*diff, 8»end, 9-change

• 1

l=8\innnary , 2“timing, 3=storage, 4=vertices, 5=pde variables, 6-settings

52
• 1

n Ivl vertices triangles unorm enorm digits

1 1 27 50 9.40E-01 3.25E-02 1.46

2 1 27 50 1.40E+00 6.81E-02 1.31

3 1 27 50 1.84E+00 7.25E-02 1.40

1 2 45 90

2 2 45 90

3 2 45 90

l=summary, 2=timing, 3=storage, 4=verbices, 6=pde variables, 6=settings

• 2

timing distribution

seconds percent

plot 0.000 0.00

assembly 7.980 50.44

factor 2.480 15.68

mg 1.210 7.65

rhs check 3.080 19.47

refine 1.070 6.76

total 15.820 100.00

l=8ummary . 2*timing, 3=storage, 4=vertices, 5=pde variables, 6=settings

• 3

45 vertices used.

1710 vertices allowed.

97 trieingles used.

4559 triangles allowed.

l=summary , 2-timing, 3=storage, 4=vertices, 5=pde variables, 6=settings

• 0

l=print, 2=plot, 3=iter, 4=refine, 5=corap, 6=save, 7=diff , 8=end, 9=chaoige

• 2

l^triamgles, 2*contours, 3=surface, 4=profiles, 5=flow lines, 6=settings

• 1

enter level

• 1

plot 0.150 sec. enter level

• 2

plot 0 . 220 sec

enter level

• 0

l=triauigles , 2=contours, 3=surface, 4*profiles, 5=flow lines, 6=settings

• 6

l=color, 2=erase, 3=magnify, 4=window, 5=file print/plot, 6=transform

• 6

enter index of transformation

• 2

transformation number 2 as PDE number 1 l=color, 2=erase, 3=magnify, 4=window, 5=file print /plot, 6=transform

• 0

53
l=triangles, 2=contours, 3=surface, 4=profiles, 5=flow lines, 6=settings

• 3

pde no., l=f\mction/2=abs(grad)

• 11

level, l=surface triangles, 2=surface contours, 3=surface grids

• 21

viewing direction — nx, ny, nz

• 1-21

min and max values -1 .334548E-03 O.OOOOOOE+00

plot 0.940 sec. pde no., l=function/2=abs(grad)

• 11

level, l=surface triangles, 2=surface contours, 3=surface grids

• -2 2

number of contours

• 9

viewing direction — nx, ny, nz

• 1-21

min and max values -1 .334648E-03 O.OOOOOOE+00

contour values

-1.201E-03 -1.068E-03 -9.342E-04 -8.007E-04 -6.673E-04

-5.338E-04 -4.004E-04 -2.669E-04 -1.336E-04

plot 0.460 sec. pde no., l=function/2=abs(grad)

• 11

level, l=surface triangles, 2=surface conto\irs, 3=surface grids

• 23

number of x lines, number of y lines

• 21 11

viewing direction — nx, ny, nz

• 1-21

min and max values -1 .334648E-03 O.OOOOOOE+00

plot 1.610 sec. pde no., l*function/2=abs(grad)

• 00

l=triangles, 2=contours, 3=surface, 4*profiles, 5=flow lines, 6=settings

• 0

Imprint, 2=plot, 3=iter, 4=refine, 6=comp, 6=8ave, 7=diff, 8=end, 9»change

• 8

54
2

Figure 1

Figure 2
Figure 3
Figure 4

c
Figure 5


NBS-n4A (REV. 2 "ac) U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET fSee instructions)

1, PUBLICATION OR

REPORT NO.

NBSIR 86-3411

2. Performing Organ. Report NoJ 3. ?'_olication Date

rULY 1986

4. TITLE AND SUBTITLE

B2DE - A Program for Solving Systems of Partial Differential Equations

in Two Dimensions

5. AUTHOR(S)

James L. Blue

6. PERFORMING ORGANIZATION (If joint or other than N BS, see in struction s)

national bureau of standards

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

7. Contract/Grant No.

8. Type of Report & Period Covered

9.

SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS ^Street. City, State, ZIP)

10.

SUPPLEMENTARY NOTES

i~ ! Document describes a computer program; SF-185, FIPS Software Summary, is attached.

11.

ABSTRACT (A 200-word or less factual summary of most significant information. If document induces a significant

bibliography or literature survey, mention it here)

B2DE is a program for solving systems of nonlinear elliptic partial differential

equations (PDEs) in two dimensions. The program is a collection’ of modules

with an interactive driver. Many types of interactive graphics plots are included

.

Users may modify the driver, and may be able to construct a "black box" program

for a restricted class of PDEs.

12. B2DE is available from the author at the above address.

KEY WORDS (Six to twelve entries; alphaoetical order; capitalize only proper names; and separate key words by semicolon si Adaptive mesh; elliptic; graphics; interactive; nonlinear equations;

partial differential equations

13. AVAILABILITY 14. NO. OF

PRINTED PAGES

, X] Unl imited

! i For Official Distribution. Do Not Release to NTIS 62

1

i Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.

20402. 15. Price

Y Order From National Technical Information Service (NTIS), Springfield, VA. 22161 $9.95

JSCOVM-OC 60O-O.0
ttffl33i-s:e soJtvAoa Itii 4

- SflSS -|

mmm ;^§U . «nolac»aia OWX Hi ---'

• 1

-1 . - ^

,***«>fr >i# TM'a 4AI i^TAO 3WU^J»«k"

ijt tfwi yn »*«S ) "t 3

ajritivfltoH^ ajtr >

.J'csksmt

sj , to«<k0M»n ^AH*^h JA*«(TA« ;r

Wt*! -gC

3S o>^

» 1

P^«0^r.3^« «*t‘!'2J?i^if|

iirp iSitivoi^ 5o i30i:>A2>JUo£i a ti •(!? , »ooiintiftlk

^

^n»>^ .>: raft

«i3 •!<*» •? :•« ^c^i ,'»»vla^ »^? txxhjom Xf

r_ J/ ‘ 4TM* _ OAtf' % ^ . . 1 f > A V <*k YJKft Y \ Jt . ^4> * # • i o« *

. mildv* ptixi* Alii JiA «Ka ir ^ *ldalUv# 1(

Jj ' ..rj !

..V'' V«v 4^ :>»%*€ •> * - • * ^ 9 *^‘

' t'' . . - - -

' ""

’

' '- ** I " ***^!—^!*8>Bd|

\-»U*

tItwrOoli jZ. ,

I . ' » « mC .1^ *•*»' #«C'*-« Vff^ »f A

t ,« t Cl ^V ***r i A' »•# »»•/ T H>

1 1 II liwi
)