Acoustic FDTD Solver

Preface
This software is the result of my experiment to create a flexible numerical simulator. It solves a class of wave problems in acoustics. It consists the major components of a numerical software: a modeller in which scatterers can be constructed by stacking solids on top of each other, a gridder to generate the rectangular mesh, a solver that implements the numerical algorithm, some visualization and post processing capabilites.
Many efforts were spent on choosing the right architecture and the right tools. Objectoriented design patterns were used in the modeller. In the meantime, the numerical algorithm was implemented in a strictly procedural way. The XML format was chosen as the input file format as it was both human friendly and machine friendly. The H5 format was chosen due to the vast amount of data that could be potentially generated. I stayed with Matlab for the visualization and plotting needs.
There are many improvements that can be added to the software, such as handling inhomogeneous media, treating curve boundaries, parallelizing using OpenMP and MPI, optimizing the use of auxillary variables in the current implementation, generating animated plots, etc.. These will be considered when it comes with a research need.
In the meanwhile, if you want to comment on the software, want to see some features, or encounter any issue, feel free to inform me.
Thank you for the interest in this work.
(g.kev) in.z hu@ie ee.o rg
Release History
If time permits, I am flexible to provide bug fixes and new features. Feel free to drop me a message. I expect to have a release every half a year.
Table of Contents
Acoustic FDTD Solver (AC2D) is a software to simulate acoustic wave propagation in lossy media in two dimensions. It uses the finitedifference timedomain (FDTD) method to solve the wave equation. It has the following features
Simulate the 2D acoustic wave propagation in lossy media. The particle velocity is in the plane of propagation, i.e. perpendicular polarization.
Implement a hard point source and a hard plane wave source at an oblique angle.
Support the pressurerelease boundary condition, the rigid boundary condition, and the periodic boundary condition. The infinite boundary condition is realized through the complex frequency shifted perfectly matched layers (CFSPML) [Gedney2005].
The software provides an executable called "ac2d", which reads a script file in the XML format as the input argument and stores the field values in an output file in the H5 format. The software also provides an easytouse C++ interface for experienced developers to explore its functionality.
The software can be used to solve the following practical problems
Solve for the wave scattered by acoustic objects due to a plane wave or an isotropically radiated wave.
Find the reflection coefficient and the transmission coefficient of stratified media.
Can be configured to run subband FDTD simulations to simulation wave propagation in linearly dispersive media [Zhao2006].
The software has been test and validated on Linux and Mac OS X. It depends on the following external software
RapidXML to parse the input script file.
HDF5 library to store data in the output file.
Boost program_options library to parse the command line arguments.
Matlab to generate the figures in the validation examples.
Doxygen to generate this documentation.
The software is placed under the open BSD license.
In linear acoustic theory [Kinsler2000], the two field quantities, the acoustic pressure denoted by and the particle velocity denoted by , are related by the continuity equation, which is given by
and the simple force equation, which is given by
where
Due to the large medium intrinsic impedance ( ), we introduce the quantities normalized with respect to the impedance , , and . The continuity equation and the simple force equation with the normalized quantities become
Converting the timedomain equations to the frequency domain (using the timeharmonic factor) gives
where the upper case letters denote the frequencydomain quantities.
The complex stretch factor in the CFSPML modifies the spatial derivatives, which is given by
where , and are the parameters in the CFSPML. The subscript can be , , or . The continuity equation, under this stretchcoordinate transformation, becomes
The simple force equation, under this transformation, becomes
One popular implementation of the CFSPML is known as the convolutional PML. Here we use the other implementation, known as the auxiliary differential equation method. Let us define an auxiliary variable (note that is spatially and temporally colocated with ), which is
and another auxiliary variable (note that is spatially and temporally colocated with ), which is
Using the newly introduced variables, the transformed continuity equation is given by
The transformed simple force equation is given by
The auxiliary differential equations that consist of the parameters of the CFSPML are given by
The corresponding timedomain equations are given by
where the lower case letters denote the timedomain quantities.
Let , and denote the spatial increments along the , and directions. Let denote the temporal increment. The acoustic pressure is located at the integer multiple of the spatial increments and the temporal increment. For instance, denotes the acoustic pressure at the location and at the time . The normalized particle velocity is located at the halfinteger multiple of the spatial increments and the temporal increment. For instance, denotes the normalized particle velocity in the direction at the location and at the time . The adopted Yee cell in the yz plane is shown below.
The centraldifference approximation of the timedomain equations are given by
The software implements the update equations in the following order: . The implementation is in ACFDTD22::execute().
Enter the following command in a terminal to install the needed libraries
$ sudo aptget install libboostprogramoptions libhdf5serialdev
The following steps are needed before compiling the software on Mac OS X
Install XCode, the integrated development environment from the App Store or from https://developer.apple.com/xcode/, free from Apple.
Inside XCode, go to XCode>Preferences>Downloads>Command Line Tools. This installs the development tools on the operating system, such as GCC/G++, make, svn, etc.
Install MacPorts, a system for developing opensource software on Mac OS X. The installer of MacPorts can be obtained from http://www.macports.org/install.php.
Open a terminal, and enter the following command to install the HDF5 library
$ sudo port install hdf518
The latest Boost source code can be obtained from http://www.boost.org/users/download/. Download the file "boost_1_51_0.tar.gz", and enter the following commands in a terminal to build Boost.
$ gunzip boost_1_51_0.tar.gz $ tar xvf boost_1_51_0.tar $ cd boost_1_51_0/ $ bash ./bootstrap # This collects the information on the operating system. $ sudo ./bjam install # This may take around 20 minutes depending on the speed of the computer.
Although this code has not been tested on Windows, it is expected to run in Cygwin, the Linuxlike environment for Windows. The HDF5 library, Boost library, and Doxygen are all available in the Package Selection Window in the Cygwin installer.
The stable release of the source code can be obtained from http://sourceforge.net/projects/ac2d/files/?source=navbar. To extract the source code from this file, enter the following commands in a terminal.
$ unzip AC2D<release number>.zip
If you have difficulty to download the file, you can access the stable release of the source code in the SVN repository hosted by SourceForge. You can see the available stable releases by entering the following command in a terminal.
$ svn ls svn://svn.code.sf.net/p/ac2d/code
You can checkout the stable release by entering the following command.
$ svn checkout svn://svn.code.sf.net/p/ac2d/code/tags/<release number> ac2d
If you want to follow the latest development of the software, you can access the code in the SVN trunk by entering the following command.
$ svn checkout svn://svn.code.sf.net/p/ac2d/code/trunk ac2d
To install the software, enter the following commands in a terminal. This generates the executable "ac2d" and places it under ac2d/bin. It also generates the library "libac2d.a" and places it under ac2d/lib.
$ cd ac2d $ make $ make install
The help message of "ac2d" is listed below.
NAME ac2d  Software to simulate acoustic wave propagation in lossy media in two dimensions. SYNOPSIS ac2d [m] [pml alpha_max=sweep xi_ratio=sweep] <file_name>.ac2d DESCRIPTION The executable reads the input script file and invokes the FDTD simulation. The output is stored in the file <file_name>.h5. OPTIONS inputfile=<file_name>.ac2d pml Estimate the optimal values of the CFSPML parameter \alpha_max and \xi_ratio. This takes the grid setting and the excitation setting (must be a GaussianModulatedSinusoid) from the input script file. It runs one simulation over a large domain and runs multiple simulations for all the combinations of the possible values of \alpha_max and \xi_ratio on a small domain. It identifies the optimal values, which yield the smallest relative error as defined in (7.135) in [Taflove2005]. alpha_max=lb:step:ub Specify the range between "lb" and "ub", in which the optimal value of \alpha_max, the shifting term, is estimated at the increment of "step". The sweep range is default to (0:0.1:1.0). xi_ratio=lb:step:ub Specify the range between "lb" and "ub", in which the optimal value of \xi_ratio, the ratio with respect to the optimal loss defined in (7.62) in [Taflove2005], is estimated at the increment of "step". The sweep range is default to (0:0.1:1.0). mem Use UPML instead of CFSPML. It sets \kappa=1.0 and \alpha=0.0 in the PML region. This reduces the memory consumption. help Display this help message. version Display the version information. AUTHOR Written by Guangran Kevin Zhu. REPORTING BUGS Report bugs to g.kevin.zhu@ieee.org COPYRIGHT Copyright (c) 20122013, Guangran Kevin Zhu
See the files placed in ac2d/validation/cpp_interface.
There are three validation examples provided in this software and they are placed under ac2d/validation. The first two examples demonstrate the use of "ac2d" and the CPP interface. The third example uses the CPP interface only.
This example demonstrates the efficiency of the CFSPML to absorb the outgoing waves. Consider a pressure pulse injected at the centre of the acoustic domain as shown in the figure below. The background has a density of 1 kg/m^3 and an acoustic speed of 1 m/s. A pressure pulse is imposed at Location O, the centre of the domain. The domain is terminated by the CFSPML in the ydirection and the zdirection. The cell size is 0.1 m in both directions, which corresponds to 10 points per wavelength at the centre frequency. Location A is 5 cells away from the PML boundary, and Location B is 5 cells away from the neighbouring two boundaries. The relative CFL number is 0.99, which translates to a temporal increment of 7.0004e02 s.
Three simulations are executed. One simulation is run with 5 PML layers on each boundary. The simulation domain has 130x130 cells. Another simulation is run with 10 PML layers on each boundary. The simulation domain has 140x140 cells. The third simulation has no PML layers. The simulation domain has 1344x1344 cells.
The pressure fields at Location A and B are recorded in the three simulations. The difference in the pressure fields between the simulations with the PML layers and the simulation without any PML layer is a measure of the ability of the CFSPML to absorb outgoing waves. This leads to the definition of a relative error given by (7.135)[Taflove2005]
The relative errors are shown in the figure below.
To run this validation example using "ac2d", enter the following commands in a terminal.
$ cd ac2d/validation/xml_interface $ bash pml_benchmark.sh
To run this validation example through the CPP interface, enter the following commands in a terminal.
$ cd ac2d/validation/cpp_interface $ make pml_benchmark $ bash pml_benchmark.sh
Inside Matlab, enter the following commands to generate the figures.
>> cd ac2d/validation >> plotCFLPMLBenchmark('xml_interface'); >> plotCFLPMLBenchmark('cpp_interface');
The two figure should be identical.
Consider a plane wave in a homogeneous medium normally incident on a cylinder as shown below. The homogeneous medium is characterized by an intrinsic impedance . The cylinder is characterized by an intrinsic impedance and a radius .
The incident field and the scattered field at the location are given by
where
The analytical expression of the field scattered by a dielectric cylinder can be found in [Balanis1989].
Consider the numerical example shown in the following figure. The background has a density of 1 kg/m^3 and an acoustic speed of 1 m/s. The cylinder has a density of 1 kg/m^3 and an acoustic speed of 0.8 m/s. It has a radius of 1 m and is centered at (9, 20). The incident pulse is a Gaussianmodulated sinusoid centered at 1 Hz. The standard deviation of the envelope is 0.5 s. The simulation domain is terminated by a 10layer PML in the ydirection, and by the periodic boundaries in the zdirection. The domain has 480x1200 cells. The cell size is 0.05 m in both directions, which corresponds to 30 points per wavelength at the centre frequency. The relative CFL number is 0.99, which translates to a temporal increment of 2.3334e02 s.
Two FDTD simulations, one with the cylinder and one without the cylinder, are executed. The difference between the fields of the two simulations gives the scattered fields at Location A and B, which are 2 m away from the centre of the cylinder. Some Matlab code is borrowed from Cylinder Scattering [CS2011], which calculates the scattered field using the analytical expression. The comparisons between the analytic solutions and the numerical solutions of the scattered fields at Location A and B are shown below. The differences between the curves are due to numerical dispersion in the FDTD method.
To run this example using "ac2d", enter the following commands in a terminal.
$ cd ac2d/validation/xml_interface $ bash plane_wave_cylinder.sh
To run this example through the CPP interface, enter the following commands in a terminal.
$ cd ac2d/validation/cpp_interface $ make plane_wave_cylinder $ bash plane_wave_cylinder.sh
Inside Matlab, enter the following commands to generate the figures.
>> cd ac2d/validation >> generateCylinderScatteredField >> plotCylinderScatteredField('xml_interface'); >> plotCylinderScatteredField('cpp_interface');
The two figures should be identical.
This example demonstrates the obliquely incident plane wave reflected by a half space. Consider a plane wave in a homogeneous medium characterized by an intrinsic impedance , obliquely incident at an angle on a half space characterized by an intrinsic impedance . The plane wave is reflected at an angle and refracted at an angle as shown below.
The particle velocities are in the plane of incidence. This corresponds to the perpendicular polarization in electromagnetics. The reflection coefficient, denoted by , is given by
and the transmission coefficient, denoted by , is given by
The angle at which the reflection coefficient is 0 is known as the intromission angle. This angle is the same as the Brewster angle of no reflection in electromagnetics. The Brewster angle for the perpendicular polarization in electromagnetics is not easily observed, since naturally occurring materials usually have very similar permeabilities. In acoustics, if the density and the velocity of the background and the half space satisfy the following relationship [DAguanno2012]
then the acoustic intromission angle is given by
Consider the numerical example shown in the following figure. The background has a density of 1 kg/m^3 and an acoustic speed of 1 m/s. The medium has a density of 1.75 kg/m^3 and an acoustic speed of 0.81228 m/s. A timeharmonic plane wave with a unity amplitude at 1 Hz incidents at an angle, and the perfectly matched layer are used at the four boundaries to terminate the simulation domain. There are 30 points allocated per wavelength. The CFL number is 0.99, which translates to the temporal increment of 3.500e02 s.
At each incident angle, two FDTD simulations, one without the acoustic half space and one with the acoustic half space, are executed. The difference in the field values recorded along the line sensor in the two simulations is the refracted wave. The reflection coefficient is the magnitude ratio between the reflected wave and the incident wave. The transmission coefficient is the magnitude ratio between the refracted wave and the incident wave. These angledependent coefficients are shown below. It is observed that, when the incident angle is 60 degrees, the reflection coefficient is 0.
To run this validation example, enter the following commands in a terminal.
$ cd ac2d/validation/cpp_interface $ make plane_wave_reflection $ bash plane_wave_reflection.sh
Inside Matlab, enter the following commands to generate the figures.
>> cd ac2d/valication >> plotHalfspaceTransmissionCoefficient('cpp_interface');
Currently, this example is not supported by "ac2d".
This software document in the HTML format can be generated by entering the following commands in a terminal.
$ cd ac2d $ make doc $ firefox doc/html/index.html
[Balanis1989] C. A. Balanis, Advanced Engineering Electromagnetics. New York, NY: John Wiley & Sons, Inc., 1989.
[CS2011] G. K. Zhu. Cylinder scattering. Online. [Online]. Available: http://www.mathworks.com/matlabcentral/fileexchange/30162cylinderscattering
[DAguanno2012] G. D'Aguanno, K. Q. Le, R. Trimm, A. Alù, N. Mattiucci, A. D. Mathias, N. Aközbek, and M. J. Bloemer, "Broadband metamaterial for nonresonant matching of acoustic waves," Scientific Reports 2, no. 340, pp. 1–5, 2012.
[Gedney2005] S. D. Gedney, "Perfectly matched layer absorbing boundary conditions," in Computational Electrodynamics: the FiniteDifference TimeDomain Method, A. Taflove and S. C. Hagness, Eds., 3rd ed. Boston, MA: Artech House, 2005, ch. 9.
[Kinsler200] L. E. Kinsler, A. R. Frey, A. B. Coppens, and J. V. Sanders, Fundamentals of Acoustics. New York, NY: John Wiley & Sons, Inc., 2000.
[Zhao2006] Y. Zhao, Y. Hao, A. Alomainy, and C. Parini, "UWB onbody radio channel modeling using ray theory and subband FDTD method," IEEE Trans. Microw. Theory Tech., vol. 54, no. 4, pp. 1827–1835, Apr. 2006.