Sven Ivansson Feb. 27, 1998
National Defence Research Establishment (FOA)
S-17290 Stockholm, Sweden
email: sveni@sto.foa.se
General information
-------------------
The rpress package is one of two packages for computation of seismo-acoustic
wave-fields in range-independent fluid-solid media that have been developed
at FOA. (The other is Ilkka Karasalo's XFEM.) As available here, the rpress
package contains five programs:
1) rpress for computation of the pressure field in a fluid region
in a fluid-solid medium primarily by wavenumber integration
(information file: rpress.txt in this directory rpress)
2) rpressfw for computation of the full wave-field in a fluid-solid
medium primarily by wave-number integration
(information file: rpressfw.txt in this directory rpress)
3) rmodfaut for computation of the modal wavenumbers and group velocities
using a winding-number integral method
(information file: rmodfaut.txt in this directory rpress)
4) rmodpfw for computation of mode shapes as functions of depth
including appropriate normalization or excitation factors
(information file: rmodpfw.txt in this directory rpress).
5) rmodfgr for computation of a dispersion curve, the modal wavenumber
of a particular mode is tracked as a function of frequency
(information file: rmodfgr.txt in this directory rpress)
The fluid-solid media that are handled consist of a fluid region interspaced
between two (possibly vanishing) solid regions. The lower and upper
boundaries may independently be either free, rigid, or of "homogeneous
half-space" type (either fluid or solid). Note, for example, that
"solid-fluid-solid" media with surrounding half-spaces of either fluid or
solid type can be handled. In addition, an entirely solid medium (without
fluid) can be handled by all programs except rpress itself. Hence, the
programs can be used not only for applications in underwater acoustics but
also in plate acoustics and seismology. The only sources that have been
implemented are (arrays of) symmetric point sources, but extensions to
double-couple sources etc. would of course be straight forward. (Relevant
formulas can be found in section 2 of Ref. 0.)
Although rpress and rpressfw were originally developed for wavenumber
integration, the field can actually be computed as a flexible mixture of
wavenumber-integral and normal-mode contributions. The procedure is to
start by computing the desired modal wavenumbers using rmodfaut, and to
continue by computing the corresponding mode shapes and excitation factors
using rmodpfw. Modal synthesis for the wave-field can then be performed
by rpress and/or rpressfw using rmodpfw outfiles as input.
The change of mode shape with frequency can be studied by running rmodpfw
using rmodfgr outfiles as input.
The underlying theory and the computational techniques are described in
the references. Basically, the compound-matrix method is used to solve
boundary-value problems in the vertical direction for fixed horizontal
wavenumbers (Ref. 0; Ref. 3, Ref. 6, Ref. 7, Ref. 8). Adaptive numerical
techniques with error control are used to do wavenumber integration and
search for modal wavenumbers (Ref. 0; Ref. 1, Ref. 2). Accuracy and
error control are given high priority, and efficiency has been pursued
with this constraint in mind. With the exception of fft transformation to
the time domain, the computations are performed in double precision.
The development and the usage of the programs have inspired some theoretical
work on wave propagation in range-independent fluid-solid media (Ref. 0;
Ref. 4, Ref. 5, Ref. 9).
The name rpress is short for "Reflectivity PRESSure". In fact, Rainer
Kind's reflectivity program for surface seismics was an important point of
departure for the work leading to the rpress programs. Another important
point of departure was Ilkka Karasalo's suggestion to amend Kind's program
with his general-purpose subroutine for adaptive numerical integration.
Implementation, some remarks
----------------------------
The directory rpress/slatec contains a "top-level" subroutine airy4e.f and
some SLATEC routines (zairy,zbiry) which are used to compute Airy-functions.
To speed up the computations, please read the header to airy4e.f, with
'note 1' and 'note 2', and take appropriate action as desired.
It might be appropriate to choose another set of data statements for
machine constants in the files i1mach.f, r1mach.f, d1mach.f in the
directory rpress/slatec. At the moment, machine constants for HP 730
have been chosen.
In the directory rpress/ik, there are a few C routines that are used to
achieve dynamic allocation of memory. Except for this (and possibly a few
details), standard Fortran 77 has been used. Each line may contain more
than 72 characters, however, and compilation should be done with the +es
option (or something equivalent). On a SUN workstation, I think it is
necessary to add a _ character to the C subroutine names in the ".c" files
in order to achieve successful communication between Fortran and C.
When linking each of the five programs, it is important that only the
appropriate directories are included, as specified in the corresponding
information file. (There is otherwise is risk of getting the wrong routine.
The names cdetsub and dhcomp appear as entry points in files in two
directories, rpress/cpm/dhcomp and rpress/cpm/fw.)
In accordance with the Fortran type rule, integer parameters are
characterized by an initial character among i,j,k,l,m,n. File names are
read as strings with 32 characters.
Test examples for running the programs are given below and in the other
information files. It is of course convenient to copy the relevant input
data to a file to which the terminal input is redirected. Some figures for
the test-example results that do not appear in the references have been
included in postscript files in this directory rpress. They are vspfig.ps
(see the information file rpressfw.txt) and slowfig.ps,slowfigp.ps (see the
information file rmodfgr.txt).
Error messages are provided in a somewhat rudimentary fashion. In general,
"STOP number" statements rather than explanatory text messages appear in the
Fortran code. Hence, one has to search the code for the appropriate number
to find the reason for the failure.
Subroutine directories
----------------------
The source code for the different programs and subroutines can be found in
the following directories:
rpress/cpm
rpress/cpm/single
rpress/cpm/dhcomp
rpress/cpm/fw
rpress/cpm/fw2
rpress/int
rpress/rut
rpress/bess
rpress/slatec
rpress/ik
The routine dpclin in rpress/int is not used in ordinary applications of
the programs. It is a Dormand-Prince ODE system solver provided by Leif
Abrahamsson and taken from the book "Solving Ordinary Differential
Equations: I. Non-stiff problems" by Hairer, Noersett and Wanner.
Some LAPACK routines to obtain machine constants appear in the files
dlamch.f and lsame.f, also in the directory rpress/int. They have been
provided by Ilkka Karasalo.
Bessel-function routines appear in rpress/bess. The routines in besbod.f
have been developed by Anders Bostrom (Chalmers University of Technology,
Gothemburg). The routines in besrat.f have been provided by Ilkka Karasalo,
they are taken from the book "Computer Approximations" by J.F. Hart et al.
Airy-function routines from SLATEC appear in rpress/slatec.
The routines in rpress/ik, for dynamic allocation of memory, have been
obtained from Ilkka Karasalo.
Example
-------
Section 2 of Paper IV in Ref. 0 contains a comparison of contributions from
different parts of the field, specified in the complex slowness plane, to the
total pressure field for a particular case. In particular, normal-mode and
branch-cut parts were considered. We will now show the required input to our
programs for a computation of the total field using a representation in
terms of normal-mode and hyperbolic branch-cut contributions. To simplify
the medium infile, combi.m below, a constant velocity profile is chosen in
the water. Four frequencies are considered: 50, 100, 150 and 200 Hz. Three
of our programs must be run in sequence:
(1) rmodfaut ---> (2) rmodpfw ---> (3) rpress
with each program run directly for all four frequencies. The required
infiles to a following program will be available as outfiles from the
previous program, with the appropriate frequency indicated by the extension
to the file name. For an explanation of the input parameters etc., we refer
to the information files rmodfaut.txt, rmodpfw.txt and rpress.txt,
respectively, all in this directory rpress.
At first, we provide the medium infile, which is given the name combi.m.
(In the information files, it is called "film".) Note that the source depth
is in the middle of the water column, and that there are ten receiver
depths to prepare for a 10x10 receiver grid in depth-range (cf. section 2
of Paper IV in Ref. 0).
0.000 0.100 5 10 ! matching depth = source depth = 0.050 km
0
0.000 0.000 0.000 0.0 0.0
0.000 0 0
1.438 0.000 1.000 0.0 0.0 0 0
0.100 0 0
1.438 0.000 1.000 0.0 0.0 10 10
0.100 0 0
1.460 0.834 1.300 0.30 0.68 0 0
0.115 0 0
1.460 0.834 1.300 0.30 0.68 1 0
0.115 0 0
4.000 2.309 2.620 0.36 0.81 0 0
(1) The first program to be run is rmodfaut, for determination of the modal
slownesses, and the required input is shown below.
To specify hyperbolic branch-cuts, "icuttyp0" is put to 1.
The modal slownesses in the upper half-plane are desired, and a particular
search rectangle in the complex slowness plane is specified by the
parameters "upreal1(),upreal2(),upimag1(),upimag2()". Since the maximum
medium slowness is 1/0.834 = 1.199 s/km, we do not expect modal slownesses
with real parts outside the chosen interval ( -2 s/km , +2 s/km ). Given
our four frequencies, and the particular source-receiver ranges to be
considered (see the input to the program rpress below), contributions from
normal modes with imaginary slowness part greater than the chosen 0.5 s/km
should be negligible.
combi.m ! film
1 1.0 0 ! mhq0,freqmhq0, iazimi
1 ! icuttyp0
1 ! isheettyp0
1D-5 1D-7 .1 1.05 .0 1.0D10 ! epsglobln,epsglob,epsfakt,hfaktmax,hoff,eps..
1D-7 0.25 ! epsfsub,hstepfsubpre
0.0 -2 0 2 7 ! xytol0, ixtyp0,nrus0,irus0,jrus0
50.0 200.0 4 ! freqmin,freqint,nfr
1.0 1.0 ! aladauer,preala
0 ! iurot
1 ! nrect
-2.0 2.0 0.0 0.5 ! upreal1(),upreal2(),upimag1(),upimag2()
1 ! ioption
0 ! icheck
-180.0 ! riktsortang
combimodu ! mainmodu
combimodux ! mainbox
0.0 0.05 5.0 ! ddsupmin,ddsupmax,ddsupfaktmax
0.12 0.4 0.75 ! derivdrelmax,drelmax,drelgoalfakt
150.0 ! argdifstore
0.0 ! slimleft
0.7 0.00001 0.000000001 ! reducfakt,ddelt,xytol
(2) The second program to be run is rmodpfw, with its "iuttyp=0" option,
for computation of the excited modes. The required input is shown below.
To specify hyperbolic branch-cuts, "icuttyp" is put to 1.
Since pressure output is desired, "kfwtyp" is put to 4.
combi.m ! film
1 1.0 0 ! mhq0,freqmhq0, iazimi
1 ! icuttyp
1 ! isheettyp
1D-5 1D-7 .1 1.05 .0 1.0D10 ! epsglobln,epsglob,epsfakt,hfaktmax,hoff,eps..
1D-7 0.25 ! epsfsub,hstepfsubpre
50.0 200.0 4 ! freqmin freqint nfr
1.0 1.0 ! aladauer,preala
0 ! iurot
combimodu ! mainmodu
0 ! iuttyp
4 ! kfwtyp
combimodp ! mainmodp
1 ! inolx
1 ! iystab
(3) The third and final program to be run is rpress, for computing the
branch-cut integrals and synthetizing the total field, including the
important normal-mode contributions. The required input is shown below.
To specify hyperbolic branch-cuts, "icuttyp" is put to 1.
Ten receiver ranges are specified: 5.550 km, 5.600 km, ..., 6.000 km. The
10x10 depth-range receiver grid is thereby obtained, it is the
"medium-range" grid in section 2 of Paper IV in Ref. 0.
The "ibessarr(),npuarr()" data of the two field parts ("nparts" is set to 2)
are important. As for the Bessel functions, H0(1) should be used in this
case with normal-mode and branch-cut contributions. Both branch-cuts are
included, and the integration is performed up to imaginary slowness parts
of 0.5 s/km (cf. the choice of "upimag2()" for the rmodfaut run).
0 10 200000 ! lx,li(top),icworkdim
combi.m ! film
1 1.0 0 ! mhq0,freqmhq0, iazimi
1 ! icuttyp
1D-5 1D-7 .1 1.05 .0 1.0D10 ! epsglobln,epsglob,epsfakt,hfaktmax,hoff,eps..
1D-7 0.25 ! epsfsub,hstepfsub
0.00001 2 0 2 5 ! preeps0, ixtyp0,nrus0pre,irus0,jrus0
50.0 200.0 4 ! freqmin,freqint,nfr
5.550 0.500 ! rmin,rint
combi.dat ! filo
combi.g ! filg
1 ! idb
combia ! mainam
combip ! mainph
0 ! intrest
1 ! ifilonpre
% ! bessfil
-1 ! iadintxfi
1.0 0 8 ! smaxpre,nelmaxpre,nbackmax
2 9 2 7 ! ixtyp,nrus,irus,jrus
0.0 10.0 430.0 ! argdiflim1,argdiflim0,argdiflim2
2 ! nparts
-2 0 ! ibessarr(),npuarr()
combimodp
0.0
-2 -12 ! ibessarr(),npuarr()
0.5
0.01 0.95 ! hmaxdh,hmaxrelwp
0.01 0.000001 ! epsh01q,preeps
0.00 0 1.0 ! hkvarpart,lunkvar,faktkvar
0 ! iconttyp
It is of course easy to compute the normal-mode and branch-cut
contributions separately. As expected, the latter are very small at the
ranges considered, but they are noticeable when high-accuracy computations
of this kind are made. At small ranges, they will definitively not be
negligible.
To compute the field using different representations and compare the
results is an instructive exercise and a useful test of accuracy. Input for
wavenumber integration (actually slowness integration) along the real axis
including "tail integrals" can be adapted from the examples in the
information files rpress.txt and rpressfw.txt (in this directory rpress).
References
----------
Ref. 0) S. Ivansson (1994):
"Seismo-acoustic wave-fields: Compound-matrix methods for
range-independent fluid-solid media". PhD thesis, Royal Institute
of Technology, Stockholm.
Ref. 1) S. Ivansson and I. Karasalo (1992):
"A high-order adaptive integration method for wave propagation in
range-independent fluid-solid media". J. Acoust. Soc. Am. 92, 1569-1577.
Ref. 2) S. Ivansson and I. Karasalo (1993):
"Computation of modal wavenumbers using an adaptive winding-number
integral method with error control". J. Sound Vibr. 161, 173-180.
Ref. 3) S. Ivansson (1993):
"Delta-matrix factorization for fast propagation through solid layers
in a fluid-solid medium". J. Comput. Phys. 108, 357-367.
Ref. 4) S. Ivansson (1994):
"Shear-wave induced transmission loss in a fluid-solid medium".
J. Acoust. Soc. Am. 96, 2870-2875.
CORRECTION: In equations (18) and (21), 1/4 should be changed to 1/2.
Ref. 5) S. Ivansson (1995):
"Source function to generate an individual mode in a fluid-solid
medium". J. Sound Vibr. 186, 527-534.
Ref. 6) S. Ivansson (1997):
"The compound-matrix method for multi-point boundary-value problems".
Z. angew. Math. Mech. 77, 767-776.
Ref. 7) S. Ivansson (1998):
"The compound-matrix method for multi-point boundary-value problems
depending on a parameter". Z. angew. Math. Mech. 78, 231-242.
Ref. 8) S. Ivansson (1998):
"Comment on 'Free-mode surface-wave computations' by P.W. Buchen and
R. Ben-Hador". Geophys. J. Int. 132, 725-727.
Ref. 9) S. Ivansson (1998):
"A class of low-frequency modes in laterally homogeneous fluid-solid
media". SIAM J. Appl. Math. 58, ?-? (in press).