Adaptive time stepping for generalized convolution quadrature

Through the current work, we aim to improve the boundary element method (BEM) in acoustics. The proposed adaptivity is based on a time‐domain boundary element formulation. In this setting, an elegant method for solving the acoustic problem is through the transformation of the underlying partial differential equation on to the boundary of the domain as a space–time boundary integral equation. However, accurate modeling of temporal behavior is crucial, requiring careful consideration of the time step size to ensure numerical stability and capture the dynamics of the system accurately. The boundary integral equation in the time domain has a convolution in time and the generalized convolution quadrature method (gCQ) is introduced, providing a framework for numerically evaluating the convolution integral within the boundary integral equation. The gCQ method under consideration uses a higher order time stepping method from the Runge–Kutta family. The gCQ method allows for adaptive time stepping, enabling the refinement of the time step size to focus computational resources where they are most needed. The current work shows that an approach similar to that of the ordinary differential equations can be included in the gCQ to introduce an adaptive control of the time step size in the 3D boundary element formulation. Numerical experiments are conducted to validate the effectiveness of the proposed approach. The results demonstrate improved solution resolution, particularly in capturing steep changes and localized variations through the adaptive time stepping scheme. The proposed method shows promising performance, dynamically adjusting the time step size based on the evolving solution, thereby enabling efficient and accurate computations in BEM simulations of radiation phenomena.

In the context of time domain simulations, the accurate modeling of temporal behavior is essential.This necessitates careful consideration of the time step size to ensure numerical stability and capture the dynamics of the system accurately.The work by Mansur [1] provided one of the earliest examples of time-domain BEM formulations, showcasing the potential of this method in capturing time-dependent radiation phenomena.A critical component of time-domain BEM is the numerical solution of the convolution integral within the boundary integral equation (BIE).The quadrature rule developed by Lubich in the 1980s has been used to approximate the convolution integral [2].However, this method is limited to a constant time step size.To address this limitation, the generalized convolution quadrature (gCQ) method was introduced by Lopez-Fernandez and Sauter [3].The gCQ method provides a framework for numerically evaluating the convolution integral with the flexibility to accommodate adaptive time stepping.The use of adaptive time stepping can significantly reduce the computational cost by ensuring that the computational resources are focused where they are most needed, leading to more efficient simulations.In a master thesis at the University of Rome under the supervision of Lopez-Fernandez [4], it has been shown that adaptivity can be introduced in the gCQ to solve the retarded integral equation.
Motivated by the potential of the gCQ method and the need for adaptive time stepping in time-domain BEM, this paper aims to introduce a novel approach that combines a Runge-Kutta-based gCQ within a space-time boundary element formulation.This combined framework allows for accurate and efficient simulation of radiation problems by integrating the advantages of BEM, gCQ, and adaptive time stepping.

PROBLEM SETTING
Let Ω ⊂ ℝ 3 be a bounded Lipschitz domain in the three-dimensional space ℝ 3 , characterized by its boundary denoted as Γ ∶= Ω.The boundary Γ could be further subdivided into disjoint sets where the Dirchlet and Neumann boundary conditions could be prescribed based on the problem.The current focus of our investigation deals with the scalar wave equation where (, ) is the acoustic pressure and  denotes the speed of the traveling wave.We are solving the wave equation for a fixed final time  > 0 with the boundary and initial conditions Our primary objective is to solve this wave equation by employing the BEM in time domain.In our example, we focus on a Dirchlet problem, while acknowledging that the proposed methodology is applicable to both types of boundary conditions.

Boundary element formulation
For the numerical solution of the wave equation, we use an indirect boundary element formulation with the single-layer potential as the solution ansatz given by Φ ∶ Ω × (0, ) → ℝ is an unknown density function and  * ∶= (−‖‖)

4𝜋‖𝑥‖
is the fundamental solution of the wave equation that refers to the solution of the underlying partial differential equation when a Dirac distribution is used as the inhomogeneity.Upon applying a suitable trace operator, the single-layer operator can be restricted to the boundary where the Dirchlet boundary conditions (  ) are prescribed.This gives us a BIE that can now be solved for the density function Φ.In the indirect BEM, the unknowns are not directly computed on the boundary of the problem domain.Instead, the BIE is used to represent the unknowns in terms of a density function defined on the boundary.The density function is then later used in Equation (3) to solve for the quantity (, ).
For space-time problems that do not allow analytical solutions, it is necessary to have a discrete representation of the underlying continuous problem.For the spatial discretization using BEM, the first step is to represent the boundary Γ by Γ ℎ such that Γ ℎ = ⋃  =1   (where   denotes the geometrical elements, e.g., triangular mesh as done in this work).The unknown field quantities are then approximated by a linear combination of some shape functions.For the density function in the indirect method, it is Φ ℎ = ∑  1 =1 Φ  ()  ().The spatial shape functions are then inserted in the BIE and a Galerkin or collocation method is applied to get the semidiscrete equation   =  * .
The obtained semidiscrete equation is a convolution integral that will be numerically solved by a quadrature rule called the gCQ method.In (4),  * is called the convolution kernel and the gCQ method uses an inverse Laplace transform to represent this convolution kernel.As the gCQ is derived, an associated ordinary differential equation (ODE) emerges.The temporal discretization is then based on a solution of this ODE using a time stepping scheme.Such an ODE within the gCQ framework facilitates the utilization of theories from ODEs for deriving the error estimates in the adaptive algorithm.For the solution of the ODE, a higher order Runge-Kutta scheme can be used for the time stepping.This leads to a higher order method in time as long as the spatial discretization is done with a higher order method as well.The mathematical analysis of a Runge-Kutta-based gCQ can be found in the works of Lopez-Fernandez and Sauter [5].
Following the aforementioned steps, the gCQ yields for a fully discretized system in space and time.The quadrature points (  ) and weights ( ) are determined with the help of Jacobi elliptic functions and detailed information can be referred from the literature [5,6].

ADAPTIVE TIME STEPPING
The main advantage of adaptive time stepping is that it allows to dynamically adjust the time step size during the computation, based on the behavior and characteristics of the solution.Fixed time step methods may require small step size throughout the computation to accurately capture all the relevant features, leading to high computational costs.Adaptive time stepping addresses this issue by selectively adjusting the time step size based on the solution's behavior.Here, we monitor the solution's local error at each time step and if the error exceeds a specified tolerance, the time step is reduced to capture fine-scale details or rapid changes.Conversely, if the error is below the tolerance, the time step can be increased to speed up the simulation.The choice of error indicator in our work is based on the classical step control for ODE's as described by Hairer et al. [7].The error is defined as the norm of the difference between a piecewise constant interpolating formula ( φ) and a piecewise linear interpolating formula ( φ) computed at the midpoint To avoid large variations in the step size, we impose a restriction and the proposed optimal time step is given by , 0.8 When the time step is accepted, we update the solution and proceed to the next iteration.Now we show the algorithm to solve the BIE using a Runge-Kutta-based gCQ with adaptive time stepping.The time steps are dynamically adjusted by controlling the error as described earlier.Algorithm 1 provided in this paper gives a pseudocode that roughly describes the key steps involved in the adaptive time stepping.In the pseudoalgorithm, we use the following notation:  The algorithm follows an iterative process until the end of the simulation time, updating the solution variable (the density function in our case) using the gCQ method.At each time step, the BIE is solved for the unknown density function.
The error is then estimated by comparing the solutions at the current and previous time steps, and the time step size is adjusted to maintain the desired level of accuracy.This process continues until the desired end time is reached.The adaptive algorithm within the gCQ framework offers the advantage of not recomputing the kernel at the quadrature points, minimizing computational overhead.Only the convolution needs to be recomputed within the loop.

NUMERICAL INVESTIGATION
The efficacy of the adaptive method is demonstrated through two numerical examples, highlighting its behavior and performance.The first example involves a cube to demonstrate the convergence characteristics of the method and a second example on a sphere addresses the accuracy of the numerical method by comparing it to an analytical solution.These are acoustic examples with the calculations being carried out for a wave speed of c = 1 m/s and an initial tolerance of  = 0.1 for the adaptive algorithm.For the numerical solution, a direct solver implemented with the HyENA boundary element library is used [8].The time descretization is achieved using a two-stage, third-order Radau IIA method from the implicit Runge-Kutta family.

Unit cube: Dirchlet problem
For the first example, we have taken a unit cube with Ω ∶= [−0.5, 0.5] 3 ⊂ ℝ 3 and the surface Γ.The boundary of the cube is discretized with triangular elements generated using the meshing software Cubit 2023.4.0 [9].A source point is chosen at the location  = (0.7, 0.7, 0.7)  and the loading is done with a traveling wave originating from here with a retarded time argument ( − ).The Dirchlet boundary data are then given by this wave having the form where () is the Heaviside function and the distance between the source point and field point is defined with  = ‖ − ‖.
A sketch of the cube along with a representation for the source point is shown in Figure 1.To observe the system's temporal evolution, the simulation is done over a time span of  = [0 s, 3 s].This particular time range is selected to ensure that a wave starting from the source point can travel the diagonal distance from the upper right corner to the lower left corner of the cube.By capturing the behavior of the wave over this time span, we can analyze its propagation characteristics and study its effects on the cube's interior.Using the indirect BEM described earlier, the density function is obtained and the values of (, ) are computed at some inner points within the cube.These computed solutions at the internal points are then used to evaluate the error for convergence.The convergence study incorporates an error measure that accounts for both spatial and temporal errors.The piecewise errors are measured using either or The equation given by ( 10) computes the piecewise error in space and the maximum value in time, whereas (9) computes the error in time by a midpoint rule.The convergence behavior is then studied using the order of convergence defined by where the errors are compared in the subsequent refinement levels based on a logarithmic scale with the base of 2. For the refinement levels, the spatial descretization is controlled such that the mesh length halves in each level.In the first refinement level, the mesh size of 0.5 m is taken.However, directly halving the time step size for each refinement level is not feasible due to the adaptive determination of the step size.Instead, the approach here is to make the tolerance  more stringent with the refinement, indirectly refining the temporal discretization to achieve a convergence in space and time.
The computed error values are presented in Table 1, with the first set of  calculated based on the  2 error and the second set based on the maximum error  max , as defined in Equations ( 9) and (10), respectively.
The Dirchlet data are approximated by linear shape functions, leading to a second-order method for the Laplace operator.For the time discretization, the two stage Radau IIA method exhibits a convergence order of 3. Considering both the factors, it is expected to have a convergence order of 2 that is confirmed in the Figures 2 and 3.
The comparison of the convergence order between an adaptive and constant step size method is also shown in Figure 3.It can be observed that the adaptive scheme exhibits slightly better convergence.However, for the second refinement level, the errors are similar for both the constant and adaptive time stepping schemes.Further inspection at this instance revealed that the adaptive scheme required fewer time steps to achieve similar error behavior as the constant step size, as shown in Figure 4.This highlights the potential of the adaptive scheme to offer improved computational efficiency while maintaining comparable levels of accuracy.

Unit sphere: Purely time depending
The main idea of an indirect BEM is to compute an intermediate value called the density function that is then subsequently used for the calculation of the desired quantity (acoustic pressure in our case).In this example, a scattering unit sphere is selected as illustrated in Figure 2.For such a configuration, it has been shown in the literature that an analytical solution In this context, the zeroth-order spherical harmonics are selected, resulting in a spatially constant solution on the surface Γ of the unit sphere.The time dependency for the Dirchlet data is given with () = −() 3  −(2−2) 2 . (13) The analytical solution is compared to numerical solutions obtained using both a constant time stepping scheme and an adaptive scheme.Figure 5 demonstrates that the numerical solutions are generally satisfactory; however, the adaptive scheme is better in capturing the steep changes within the solution more accurately and reproducing the intricate details and rapid variations of the solution.

CONCLUSION
In conclusion, this study focused on the solution of the wave equation using the BEM in the time domain.A higher order Runge-Kutta based gCQ method was successfully applied that opened up the possibility to have a higher order spacetime BEM formulation.The introduction of an adaptive time stepping scheme within the gCQ framework for a BEM setting was a significant contribution of this work.The adaptive time stepping scheme allowed for improved resolution of the solution, particularly in capturing steep changes and localized variations.Numerical results show that the proposed scheme works as expected and by dynamically adjusting the time step size based on the evolving solution, adaptive time stepping approach allows for efficient and accurate computations in BEM simulations.

F I G U R E 1
Unit cube with triangular meshes.

F I G U R E 2
Spatial descretization for sphere.F I G U R E 3 Error for the unit cube: Dirchlet problem.F I G U R E 4 Number of steps comparison for a constant step size and an adaptive step size schemes.F I G U R E 5 Comparison of density function for numerical and analytical solutions.existsbased on the eigenpairs of the boundary integral operator.If we represent the sphere with Ω ⊂ ℝ 3 and    denotes the spherical harmonics, then for a given right-hand side   (, ) = () 0 0 , the analytical solution for the density function of the single-layer approach is given by[10] ,   ← quadrature nodes and weights for  = 1, … ,   quadrature points 5: Δ  and Δ  ← determine the maximum and minimum step size  −1 (  ) = ( − Δ −1   ) −1 ((   −1 ⋅  −2 (  )) + Δ −1  −1 ) ← Perform the Runge-Kutta step • Δ for the time step sizeThe following parameters are prescribed by the user:• Δ max → maximum allowed time step size • Δ min → minimum allowed time step size •  → the prescribed tolerance below which the error should be kept •  → maximum number of rejected time steps at each iteration Error estimates and convergence order.
TA B L E 1