Sven Ivansson Oct. 15, 2022 (originally Feb. 27, 1998) Swedish Defence Research Agency SE-16490 Stockholm, Sweden email: sven.ivansson@gmail.com 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 the Swedish Defence Research Agency. (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) NOTE that, in addition to the original double-precision code, * there is a quadruple-precision version * there is a quadruple-precision version for SH waves called rmodfautSHf (see the end of rmodfaut.txt) 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). A more general presentation of underwater sound propagation modeling, with many more references, appears in Ref. 10. 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 ---------------------------- To compile and link the programs, see the file INSTALL.txt The directory rpress/slatec contains "top-level" subroutines airy4e.f and airy2e.f, and some ACM routines (aiz,biz) which are used to compute Airy- functions. To speed up the computations, please read the headers to airy4e.f and airy2e.f, with 'note 1' and 'note 2', and take appropriate action as desired. Except possibly for a few details, standard Fortran 90 has been used. Each line may contain more than 72 characters, and compilation should be done with this in mind. 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, often 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/linu and also, concerning quadruple precision, in rpress/cpmQ rpress/cpmQ/singleQ rpress/cpmQ/dhcompQ rpress/intQ rpress/rutQ rpress/slatecQ 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 ACM appear in rpress/slatec. 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 Since one "npuarr()" vanishes in this example, an input file UPVAL.DAT is needed (see the information file rpress.txt). To include all available modes in the computations, UPVAL.DAT can consist of the single line 0 ! iupval 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, 1462-1508. Ref. 10) S. Ivansson (2017): "Sound propagation modeling", in Applied Underwater Acoustics, chap. 3, 185-272, eds. L. Bjorno, T. Neighbors, and D. Bradley, Elsevier.