NetNews Usenet Archive 1992 #16/NN_1992_16.iso/spool/comp/lang/fortran/2823

Internet Message Format | 1992-07-24 | 9.1 KB

Path: sparky!uunet!olivea!decwrl!sun-barr!cs.utexas.edu!torn!cunews!nrcnet0!cu23.crl.aecl.ca!wl.aecl.ca!harrisp
From: [e-mail]
Newsgroups: comp.lang.fortran
Subject: SUMMARY: Partial Diff. Eqn's. Thanks to all
Message-ID: <[e-mail]>
Date: 24 Jul 92 18:39:02 GMT
Sender: [e-mail] (USENET News System)
Organization: AECL RESEARCH
Lines: 209
Nntp-Posting-Host: wc4.wl.aecl.ca

From: [e-mail]   (Phil Harris)
Date: 24-JUL-1992
RE:   Solving partial differential equations numerically (SUMMARY)

 I originally posted in article <[e-mail]>

>From: [e-mail]   (Phil Harris)
>Date: 22-JUL-1992
>RE:   Solving partial differential equations numerically
>
>Can anyone recommend a fortran callable  package for solving pde's and systems
>of pde's. I really need the source code, preferably written as standard
>Fortran 77 as I want this to run on a PC (80486) using Lahey Fortran. Some
>documentation would be nice also. I'm looking for something that is fast,
>since the model I am developing has to run for several months of real time.
>I wrote a simple program to solve the equations using explicit
>finite differences which works just fine, but I need to use very small time
>steps for it to remain numerically stable, and takes almost a minute of
>execution time on my '486 to simulate the results of a minute of real time.
>I am currently working on a code using implicit (Crank-Nicholson) finite
>differencing, but am having some trouble getting it to execute correctly.
>I'm hoping that once I get this working, I will be able to use larger time
>steps.
>
>The equations I need to solve are fairly simple:
>
>dQ(z,t)
>-------  = a1*( Const1(z) - Q(z,t) )*C(z,t) - a2*Q(z,t)
>  dt
>
>dC(z,t)                dQ(z,t)    dC(z,t)
>-------- = -Const2 * ( ------- +  ------- )
>  dt                     dt         dz
>
>a1, a2, Const2  are independant of z and t
>Const1 independant of t, but dependant on z.
>
>Phil Harris                |  Imagination is more
>Whiteshell Laboratories        |  important than knowledge.
>Atomic Energy of Canada        |           Albert Einstein
>[e-mail]            |


I want to thank all of those who kindly responded, and who have opened my
eyes to a whole world of useful information in the NETLIB libraries. Also
thanks to those who provided several useful suggestions about implicit
methods and the application of Crank-Nicholson. I didn't mention that I
would provide a summary in my original post, but thought it might be useful
especially to those who are new to using bulletin and anonymous FTP.

            SUMMARY

Venkatesh Gopinath (Dept of Electrical Engineering,Michigan State University)
has used a set of routines from the NAG library for solving systems of PDE's
and found that they are easy to use and work well. Source code is provided,
but the NAG libraries are not free. I am  looking into their availability
for the PC.

Carl Gooch ([e-mail]) obvious knows his way around
Crank-Nicholson and implicit codes for PDE's. His advice assured me that
what I was trying should work, although he included:

>I would expect Crank-Nicholson to be fairly easy to implement, even
>with the presence of the non-linear Q*C term, which will end up
>looking something like: Q(z,t)*C(z,t) + 0.5*(Q(z,t)*dC(z,t)/dt * dt +
>dQ(z,t)/dt *C(z,t) * dt) + O(dt**2).  Two problems to keep in mind:
>first, Crank-Nicholson is only neutrally stable; numerical errors do
>not damp out.  The easiest way to fix this is to use a fraction
>slightly greater than 1/2 in expressions like the above.  That makes
>the method slightly more implicit, awhich improves the stability
>characteristics of the method at the price of dropping the time
>accuracy back to first order.

He was right in both counts, easy to implement and required lambda
not equal to 1/2. It seems my problem with my implicit code (which I
tried debugging for two weeks) was a compiler bug and not the code
itself. Might post something about this later, once I've done more
testing. Don't want to flame the supplier if I'm wrong. The implicit
code has been tested using time steps two orders magnitude larger than
those used in the explicit code, and executes about 60 times faster.
There are still some areas in the code which can be optimized, and I may
be able to go to still larger time steps without a significant loss in
accuracy.

Thanks Carl for the advice and assurance that it should work.


        NETLIB

Lot's of good advice re: NETLIB

From [e-mail] I received
copies of a number of other postings related to NETLIB, and a brief USER's
GUIDE to NETLIB. This made it obvious to me that there was a lot of useful
source code sitting out there waiting to be used.

The NETLIB routines may be obtained in two ways:

>If you have not tried netlib, try sending email to
>        [e-mail]
>with message
>        send index

or:
    FTP to research.att.com
    Password: NETLIB

I took the second route and was blown away by the amount of stuff there.


Stan [e-mail] pointed out the size of Netlib (actually
warned) and how to get access to it.

C. Wolfe ([e-mail]) also sent info on getting
access to NETLIB>

Both [e-mail] and [e-mail] sent info about
NETLIB, and a particular routing called pdecol.f.

>Use the pdecol.f Program. It uses the method of the orthogonale collocation
>. It is very usefull for great time steps. It is in the NETLIB routines.

>Algorithm 540 in TOMS is PDECOL, software for evolution of PDE's
>in time.  It works OK (but is not terribly fast).

seemingly some disagreement about speed, but both obviously found it
useful. TOMS stands for Transactions of the Mathematical Society (I believe)
and algorithms related to published articles are available on NETLIB

> I suggest you try sending e-mail to [e-mail] with the
>messages (not subject lines):

>    send index
>    send index from toms
>    send 540 from toms

I have downloaded pdecol from NETLIB, but have not yet had the opportunity
to get it up and running.

Thanks to all for the help re NETLIB. All I can say after seeing what's
there is

    So Many Interesting Things to See and Do,
    Too many deadlines to meet.


            END of NETLIB Stuff


[e-mail] pointed out that "Numerical Recipes" discussed
finite differcing schemes quite extensively, and provided some advice on
getting routines to solve the resultant equations from the finite
differencing scheme. Also, mentioned a "set of black - box routines
called SLDGL" that he has worked with and suggested they might be useful
in solving this system.


>I've also worked with a set of black - box routines called SLDGL, which
>were developed at the Computer Centre of the University of Karlsruhe,
>Germany; they could do what you want.  To get a line to the people
>there, mail [e-mail]. They also developed FIDISOL,
>a package optimized for supercomputers.

Thanks much for the advice. I had already looked at Numerical Recipes,
which helped greatly in setting up the equations. Have downloaded MADPACK
from NETLIB, but no time yet to implement. Am going to archie search for
SLDGL.

[e-mail] expressed the greatest interest in the physics of
the problem, and what the real world applications are. The problem deals
with the transport of an adsorbable ( probably not found in any dictionary
but what the hell) species through a granulated bed with constant carrier
gas flow rate. The bed contains porous particles, which have a limited
number of adsorption sites. The first equation describes the rate of change
of the concentration of the adsorbed species (Q), in which the adsorption
proces obeys a Langmuir isotherm. The constant with a z dependance is the
number  of free sites per unit length of bed. In most cases, the z dependance
is nil. However, I have to consider the case where the front parts of the
bed (small z) have been poisoned by a material which decreases the number
of active sites. The problem has been simplified by ignoring diffusion
effects.

For sufficiently long bed lengths, and in the case of a uniform distribution
of the number of active sites (ie no z dependance, no poisoning) the system
of equations will result in constant pattern conditions. This gives rise to
a wave equation for the gas phase concentration (C) and the propagation of
a shock wave solution.

An analytical solution to the system can be obtained under the conditions
that there is no z dependance in the first constant. I have not tried yet
to reproduce this solution, but it can be found in.

Thomas H. C., Ann. N.Y. Acad. Sci., Volume 49, pg 161 (1948)

            END of SUMMARY

Once again, thanks to all who took the time to respond. This was the first
time I had posted a question on the boards, and I was astounded and delighted
with the response. If I missed giving credit to anyone who responded, my
sincerest apologies.

Sorry for the length of this summary. First time I posted, but I had read
that it is good netiquette to provide a summary if the responses were many,
and may be of general interest.

Phil Harris                |  Imagination is more
Whiteshell Laboratories                |  important than knowledge.
Atomic Energy of Canada                |           Albert Einstein
[e-mail]            |