OpenPSTD v1.0 overview
What can you model with openPSTD?
openPSTD provides a visual user interface that enables to compute sound propagation in a 2D geometry. With openPSTD, the 2D geometry, bounded by edges that (partly) absorb the sound waves, can be drawn in Blender. Sources of sound and receiver positions can freely be chosen within the geometry. Sound propagation is computed in openPSTD by numerically solving the linearized Euler equations, i.e. the partial differential equations that describe sound propagation in fluid flow in the time domain. The outcome of openPSTD is a visualisation of sound propagation, by a 2D animation and by the impulse responses at the position(s) of the receiver(s). As openPSTD is an open source software, we highly encourage contributions to the code.
- openPSTD v1.0 has been launched in september 2014.
- openPSTD v1.1 has been launched in july 2015. This version is the GPU accelerated version of openPSTD v1.0
The openPSTD interface is built as a add-on for Blender, a free and open source 3D modelling application. As such, it inherits the a lot of the functionality from Blender . For a comprehensive step-by-step guide to model a geometry in Blender please see Getting started with openPSTD. For detailed instructions on how to use Blender, users are encouraged to find some of the tutorials available .
It is also possible to model a geometry without using the visual Blender interface. A guide can be found here: Command line simulations.
Geometry, medium and boundaries
A geometry in openPSTD must be constructed out of a subsequent arrangement of rectangles (subdomains). The interface, at which the boundaries of two adjacent subdomains coincide, does not necessarily have to match the full length of both subdomain boundaries. Hence, subdomains can have multiple neighbors at a domain boundary and hence local refinement does not necessarily propagate throughout the entire topology. The density and speed of sound of the propagation medium can be defined in the openPSTD settings. A restriction of the method is the minimum dimension of the subdomains, which is equal to the Grid spacing multiplied by the Window size (see openPSTD v1.0 settings). Once the geometry has been constructed, the acoustic absorption coefficient α at the edges needs to be defined. The defined absorption coefficient corresponds to a real valued normalized surface impedance according to
The absorption coefficients are frequency independent. The default setting in openPSTD is that boundary reflections are independent on the angle of incidence. A locally reacting boundary can also be chosen, see openPSTD v1.0 settings. The reflection coefficient then becomes, in approximation (see van Renterghem et al. 2013).
with θ the angle of sound incidence relative to the normal to the interface.
Source and receivers
The location of source(s) and receiver(s) can be chosen freely in the constructed geometry, see Getting started with openPSTD v1.0. Source(s) are modelled by initial pressures distributions with the following spatial Gaussian shape:
with , Δ the grid spacing in metres, c the adiabatic speed of sound in m/s, x=[x,z] the 2D coordinate vector and xs the source position. The centre(s) of the Gaussian shape(s) thus corresponds to the defined centre(s) of the sources. Although receiver positions might not correspond to a dicrete position in the mesh, spectral interpolation allows to compute the response of the source(s) in the constructed geometry at an arbitrary location in the space.
In openPSTD v1.0, the geometry is discretised by an orthogonal equidistant grid, with Grid spacing in metres equal to the spacing of the grid points. Pressure grid positions are staggered by a half grid spacing with acoustic velocity grid positions, and the velocity position are located at the edges of a subdomain: horizontal velocity components are located at a vertical edge and vertical velocity components are located at an horizontal edge. When a geometry is drawn, a grid is created. If a dimension of a drawn subdomain does not correspond to an integer multiple of the Grid spacing, the subdomain is extended such that its size (in the respective direction) becomes an integer multiple of the Grid spacing.
Wave propagation in the openPSTD software is implemented according to the extended Fourier PSTD method as presented in Hornikx et al. 2010. With this numerical method, the linearized Euler equations (LEE) are solved, which allow a detailed computation of sound propagation in moving inhomogeneous media and return the linear acoustic equations for a non-moving homogeneous propagation medium. In PSTD, the spatial derivatives arising in the LEE are computed in the wavenumber domain. This implies that a spatial Fourier transform is applied to the acoustic variables, and after an operation in the wavenumber domain, an inverse Fourier transform returns the derivatives. This way of computing spatial derivatives leads to an accurate solution requiring only 2 spatial points per wavelength. In openPSTD v1.0, the total domain is decomposed into subdomains. The advantage of this decomposition is a higher geometrical flexibility, a lower computational complexity, it promotes object orienting programming and a parallel computational architecture. The computation of the horizontal spatial derivatives in a subdomain at the pressure positions is schematically depicted in the accompanied Figure. Three subdomains are needed to compute these derivatives. Prior to computing the derivatives, the variable values are multiplied by a window function. The accuracy of this technique is reported in Hornikx et al. 2012, and increases with an increasing length of the decay of the single sided windows. Behind every edge, a boundary subdomain is included in openPSTD v1.0. This subdomain is not visible in the graphical interface in Blender. The waves propagate in the boundary subdomains and are absorbed by a perfectly matched layer (PML), see Hornikx et al. 2010.
The code is written in Python. Python takes as input a description of the boundary geometry and source and receiver positions, as it is compiled by the user interface in Blender. Python in turn outputs the resulting pressure values for the computational grid as well as at specified receiver positions, which are read back into the user interface to be visualized there.
Computed results can be accessed in two ways. First, the 2D pressure field can be viewed as a function of time. This offers the possibility to interpret sound propagation in a qualitative way. The second possibility is to export the computed impulse responses at the receiver positions for further post processing. The broad band impulse responses can then be processed to frequency domain results and octave band values can be obtained.