Based on the nodal DG-FEM Matlab implementation provided by Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications by J. S. Hesthaven and T. Warburton (2007).
PML formulation from A stable decoupled perfectly matched layer for the 3D wave equation using the nodal discontinuous Galerkin method S. J. Feriani, M. Cosnefroy, A. P. Engsig-Karup, T. Warburton, F. Pind, C.-H. Jeong (2024)
Meshes can be generated with Gmsh 1 and examples are present in the mesh folder.
To select the mesh for the simulation, change the directory of fileToLoad
.
The solver is initialized using StartUp2D
or StartUp3D
to build the spatial discretization according to the order of approximation chosen N
and the mesh characteristics. The functions to initialize the solver are contained in
fpeak
is the peak frequency of the initial Gaussian distribution for the pressure (see 2).
The function CheckPPW
makes sure that there are enough points per wavelength in each element of the mesh. It calculates also the upper frequency that can be considered (when half-power of the frequency response is reached).
Time is discretized using the low-storage five-stage fourth-order explit Runge-Kutta scheme 3.
Choose the final time of the simulation (FinalTime
) and the Courant–Friedrichs–Levy number (CFL
)
No acoustic energy should be present in the PML at the initial time (it create instabilities in the PML that could ppread to the domain of interest) so CheckEnergyPML
makes sure that the pressure distribution is narrow enough that the pressure is 0 in the elements of the PML.
The damping function PMLsigma
. In 2D you can choose the PML to implement: U-PML, M-PML, C-PML or the decoupled PML (default).
The while
loop get executed until the time has reach FinalTime
using the explicit low-storage forth-order Runge-Kutta method 4.
The if
statement inside the loop can be activate to check on the simulation as it runs every 25 time steps (the plotting of interemediary steps slows down the simulation).
At the end of each time step, the output is calculate if the variables have been initialized before the begining of the time loop.
By default: recP
, the extrapolated pressure at a receiver point and energyP
, the energy inside the domain of interest.
Footnotes
-
C. Geuzaine and J.-F. Remacle. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering 79(11), pp. 1309-1331, 2009. ↩
-
Shinichi Sakamoto, Phase-error analysis of high-order finite difference time domain scheme and its influence on calculation results of impulse response in closed sound field, 2007 ↩
-
M. H. Carpenter, C. A. Kennedy, Fourth-order 2N-storage Runge-Kutta schemes, Technical Memorandum (NASA-TM-109112), NASA angley Research Center, Hampton, VA (June 1994). ↩
-
M. H. Carpenter, C. A. Kennedy, Fourth-order 2N-storage Runge-Kutta schemes, Technical Memorandum (NASA-TM-109112), NASA Langley Research Center, Hampton, VA (June 1994). ↩