A brief description of the FENL program

N. A. Kampanis and E. T. Flouri

Institute of Applied and Computational Mathematics

Foundation for Research and Technology

71110 Heraklion, Crete, Greece

Introduction

The FENL program computes numerically the solution of the Helmholtz equation in an axially symmetric waveguide consisting of fluid layers overlying a rigid bottom, using the finite element technique. The above approach is used to simulate the propagation of sound waves due to a harmonic point source placed in the waveguide.

Mathematical formulation

The Helmholtz equation is written in cylindrical coordinates, irregular interfaces between the fluid layers may be considered, a Neumann boundary condition is applied at the bottom (which may as well be irregular), a homogeneous Dirichlet boundary condition is prescribed at the surface (which may also be irregular, e.g. in the case of a deterministic rough surface), a nonhomogeneous Dirichlet boundary condition (to be described in the sequel) at the inflow (left, vertical) boundary and a nonlocal (DtN) absorbing boundary condition at the outflow (right, vertical) boundary.

Numerical scheme

We use the standard Galerkin/ finite element method with P1-triangular elements with a weighted L2-inner product which takes into account the interface conditions. The nonlocal outflow boundary condition is suitably discretized and incorporated in the weak formulation. The finite element mesh is a general nonuniform mesh.

For details on the mathematical formulation, on the discretization technique and the practical aspects of the implementation cf. . For comparisons with coupled mode solutions of the Helmholtz equation cf. .

The FENL program

The FENL program is written in Fortan 77. It has several modules, i.e. separate Fortran codes, each of which implements a specific step of the finite element technique. These modules are:

• Renopo.f, which reads the triangulation data created by the mesh generator of the MODULEF library (the so called NOPO data structure) and fills with them the arrays required in the subsequent steps. This module is self-contained.
• Amat.f, defines the format (i.e. allocates space in the arrays defining the associated data structure) with which the matrices of the finite element discretization of the boundary value problem for the Helmholtz equation, and of the 1D eigenvalue problem from which we compute the discrete DtN map, will be stored. The nonzero structure of the lower triangular part of the matrices is stored. This module is self-contained.
• Eigen.f, solves the 1D, eigenvalue problem resulting from the discretization of the DtN map. The subroutines required for the optimal ordering on the 1D finite element matrix are from the SPARSEPACK (A. George, J. W. H. Liu, Computer Solution of Large Sparse Positive Definite Systems, Prentice-Hall, 1981) and are included. The NAG (NAG Fortran Library Manual, Mark 15, The Numerical Algorithms Group Ltd., 1991) subroutines required for the computation of the Hankel functions are included. The EISPACK (B. T. Smith, J. M. Boyle, B. S. Garbow, Y. Ikebe, V. C. Klema, C. B. Moler, Matrix Eigensystem Routines-EISPACK Guide, Springer, 1976) subroutine required is not included and the user should link with it.
• Assembly.f, computes the entries of the matrices from the finite element discretization of the Helmholtz equation. It is self-contained.
• Solve.f, solves the linear system occurring by the finite element discretization of the Helmholtz equation. First it transforms the initial format used to store the nonzero structure of the matrices, to the "compressed sparse row" format accepted by the QMRPACK (R. W. Freund, N. M. Nachtigal, QMR: A quasi-minimal residual method for non-Hermitian linear systems, Numer. Math. 60 (1991), 315-339) which is then used for the solution. This module is not self-contained. The QMRPACK has to be linked with this module. The software can be obtained via the web from http://www.netlib.org/linalg (we download the file qmrpack.tgz). In the directory MAKE we include the makefile for Solve.f under the name MakeSolve.
• Plot.f, creates the necessary output files to be used by the functions of the PDE Toolbox of MATLAB which will be used for plotting.

Mesh generation

For the mesh generation we use the mesh generator APNOXX of the MODULEF finite element library (M. Bernadou, P. L. George, A. Hassim, P. Joly, P. Laug, A. Perronnet, E. Saltel, D. Steer, G. Vanderborck, M. Vidrascu, MODULEF: A Modular Library of Finite Elements, INRIA, 1988). This library is available via the web from http://www-rocq.inria.fr/modulef. The Renopo.f module can read the data structure NOPO and extract or create the necessary information for the triangulation.

Input files

The FENL program requires the following files as input:

Intrf.inp, containing the physical parameters for the test case. These are:

1st line

FRQ (source frequency in Hz), ZS (source depth in m), C0 (reference sound speed in water in m/sec).

2nd line

ZBOT (maximum depth, at the outflow boundary placed at RMAX), RMAX (maximum range, where the outflow boundary is placed).

3rd line

CW (sound speed in the first layer-water), CB (sound speed in the second layer-bottom).

4th line

RHO1 (density of the first layer-water), RHO2 (density of the second layer-bottom).

5th line

ZSTR (the depth where the interface between the two layers intersects the outflow boundary placed at RMAX).

Files.inp, which contains:

1st line

NOPOFILE, a CHARACTER*(*) variable containing the name of the file where the data structure NOPO is stored.

2nd line

ISF, an integer variable. If set equal to 1 the nonhomogeneous Dirichlet boundary condition at the inflow boundary (left, vertical) is given in a pointwise manner, i.e. the values of the solution are prescribed at certain points on the inflow boundary, by the user. For example a normal mode code can be used to provide these pointwise values (cf. the test cases in  and ). If set equal to 0 the solution on the inflow boundary is defined analytically (one can choose the Gauss or Greene initial field).

3rd line

SFLDFILE, a CHARACTER*(*) variable, which if ISF=1 contains the name of the file where the solution is given in a pointwise manner (in pairs of the form: Z, value of solution at Z).

Technical details

From the FENL.tar.gz, after the unzip and untar procedures we have the FENL directory which contains the following subdirectories:

DATA, which contains the intrf.inp and files.inp files and the directory JCA_7_99_TC1 where the data to run Test Case 1 of  are found. These files are the QQ.NOPO, which contains the NOPO data structure and is created by the APNOXX driver of the MODULEF library, and fld.str, which contains the pointwise values defining the nonhomogeneous Dirichlet boundary condition for the solution at the inflow boundary. In this case they were produced by the COUPLE coupled mode code (R. B. Evans, COUPLE: A user’s manual, NORDA TN-332, 1986).

EXEC, where the executables must be placed (empty at first).

MAKE, where the makefile for Solve.f is found under the name MakeSolve.

QMRL, where the QMRPACK must be placed and compiled (empty at first).

SRCS, where the source files Renopo.f, Amat.f, Eigen.f, Assembly.f, Solve.f and Plot.f are found.

RSLT, where the output from each module will be directed (empty at first). The files, which will be used for plotting, are fem.out and fld.out. The fem.out file contains the neighboring nodes (in global numbering) for each elemnt. The fld.out file contains the computed value of the solution at each node of the triangulation. Both files are required by the plotting routines of the PDE Toolbox of MATLAB (cf. Partial Differential Equation Toolbox User’s Guide (for use with MATLAB), The MathWorks Inc., 1996).

Matl, where the fem.out and fld.out files from RSLT must be copied for plotting with MATLAB. The names of the new files are fem and fld, respectively, since MATLAB requires file names without a separating dot. In the subdirectory JCA_7_99_TC1 the fem and fld files for Test Case 1 of  are given accompanied by the associated MATLAB plotting program FE.m. Further, the reference solution by COUPLE is provided in the cple file as well as the associated MATLAB plotting program CM.m.

Nopo, where in the subdirectory JCA_7_99_TC1, the input file for APNOXX for Test Case 1 of  is provided.

To run a certain test case with FENL, we need first to create a triangulation of the computational domain using APNOXX driver from MODULEF library. For this we need to prepare the input file for APNOXX (find a sample in Nopo/JCA_7_99_TC1). Then we prepare the intrf.inp and files.inp (find a sample in DATA/JCA_7_99_TC1). We run the modules of FENL in the following order: Renopo, Amat, Eigen, Assembly, Solve and Plot. Then we copy from RSLT to Matl, the fem.out and fld.out files under the names fem and fld, respectively (find a sample in Matl/JCA_7_99_TC1) and plot them with the FE.m program (appropriately modified for the case at hands).

References

1. N. A. Kampanis, V. A. Dougalis, A finite element code for the numerical solution of the Helmholtz equation in axially symmetric waveguides with interfaces, J. Comp. Acoustics 7 (1999), 83-110.
2. V. A. Dougalis, N. A. Kampanis, M. I. Taroudakis, Comparison of finite element and coupled mode solutions of the Helmholtz equation in underwater acoustics, Proceedings of the 4th European Conference on Underwater Acoustics, A. Alippi, G. B. Cannelli (Eds.), CNR-IDAC, 1998, Vol. II, 649-654.