SIAM News Blog
Research

A Fast Finite Element Solver for Focused Ultrasound Applications

By

From chronic conditions like malignant tumors to movement disorders such as Parkinson’s disease, the global population experiences a diverse spectrum of health concerns. Despite significant medical advancements, current treatments are not without limitations. Therefore, the development of new treatment methods is essential to effective health management.

<strong>Figure 1.</strong> Setup of a focused ultrasound treatment, where \(D\), \(f\), and \(\Omega\) respectively denote the transducer diameter, transducer focal length, and computational domain. The purple region represents the computational domain, the pink region represents the organ, and the brown region represents the tumor. Figure courtesy of the author.
Figure 1. Setup of a focused ultrasound treatment, where \(D\), \(f\), and \(\Omega\) respectively denote the transducer diameter, transducer focal length, and computational domain. The purple region represents the computational domain, the pink region represents the organ, and the brown region represents the tumor. Figure courtesy of the author.

Focused ultrasound (FUS) is an emerging treatment modality that uses ultrasound waves to target particular tissue within the body. It has proven to be effective for cancer therapy, the treatment of brain disorders, and the delivery of therapeutic agents across the blood-brain barrier [2, 4, 6]. FUS is noninvasive—which reduces the risk of infection—and localized, so that the treatment of target tissue does not affect the surrounding tissue. Figure 1 depicts the setup of FUS treatment.

Despite the potential applications of FUS therapy, its practical use remains limited due to the difficulty of patient-specific treatment planning. For instance, healthcare practitioners must model the effect of the generated ultrasound waves on the patient, generally taking a computational approach due to the mathematical model’s complexity. Moreover, the computational model must be solvable within a matter of minutes or even seconds to prove effective for treatment planning [5]; such a feat is challenging because the inherently large domain-to-wavelength ratio leads to a high computational demand. This issue is exacerbated for high-intensity FUS applications, wherein the mathematical model’s nonlinearity leads to higher-frequency waves.

To address these computational complexities, we have developed a finite element solver for FUS applications. We consider the linear wave equation with attenuation as the mathematical model:

\[\frac{1}{\rho_0}\nabla^2p-\frac{1}{\rho_0c^2_0}\frac{\partial^2p}{\partial t^2} + \frac{\delta}{\rho_0c^2_0}\frac{\partial p}{\partial t}=0. \tag1\]

Here, \(p\) is the pressure field, \(\rho_0\) is the density, \(c_0\) is the speed of sound, and \(\delta\) is the diffusivity of sound. This equation effectively models low-intensity FUS applications that are typical in transcranial FUS simulation [1].

Next, we discretize equation \((1)\) via two methods: (i) the finite element method for spatial discretization and (ii) the explicit Runge-Kutta method for temporal discretization. For the finite element discretization, we utilize a mass-lumped scheme to yield a diagonal mass matrix that further reduces the computation and storage requirements, as only the diagonal entries of the matrix must be stored. We perform the mass-lumped scheme in the context of the nodal quadrature approach, for which the numerical quadrature points serve as interpolation points of the basis function [3]. When combined with an explicit time-stepping method, the mass-lumped scheme allows us to bypass the step of solving the linear system — consequently reducing the overall computational cost.

<strong>Figure 2.</strong> The normalized simulation time of the solver. The blue bar represents the pure FEniCSx implementation, the orange bar represents the precomputation of geometric data, the green bar represents the precomputation of geometric data with sum factorization, and the red bar represents the previous two implementations along with the use of 32-bit floating-point numbers. Figure courtesy of the author.
Figure 2. The normalized simulation time of the solver. The blue bar represents the pure FEniCSx implementation, the orange bar represents the precomputation of geometric data, the green bar represents the precomputation of geometric data with sum factorization, and the red bar represents the previous two implementations along with the use of 32-bit floating-point numbers. Figure courtesy of the author.

We implemented the solver with C++ to provide access to the low-level functions in FEniCSx: an open source software that uses the finite element method to solve partial differential equations. FEniCSx’s built-in features enable a performant numerical implementation that includes support for quadrilateral and hexahedral meshes, a high-order finite element basis function, and parallel processing capabilities. These features allow users to obtain accurate and realistic simulations of FUS applications.

One of our main goals is to improve the solver’s performance and ultimately reduce its time-to-solution. Doing so requires identifying and addressing the solver’s main computational bottleneck, which we found to lie in the assembly of the finite element vectors. To address this slowdown, we relied on three methods: (i) the precomputation of geometric data, (ii) the sum factorization algorithm, and (iii) reduced-precision floating-point numbers.

When the geometric data is precomputed, the FUS simulation mesh is typically fixed at each time step, meaning that the precomputation and storage of geometric data lowers the overall computational cost by reducing the number of floating-point operations per time step. We also exploit the tensor-product structure of the quadrilateral and hexahedral element in order to implement the sum factorization algorithm, which reduces the number of required floating-point operations for vector assembly. Finally, the use of 32-bit floating-point precision increases the number of floating-points that can be processed simultaneously, thus reducing both vector assembly time and storage requirements. Figure 2 illustrates the speedup from each implementation.

In conclusion, we have implemented a finite element solver for FUS applications that utilizes high-order basis functions, efficient algorithms, and parallel processing capabilities. We leveraged several methods to improve the solver’s overall performance and found that its accuracy is comparable to other established solvers (see Figure 3).

<strong>Figure 3.</strong> Pressure amplitude comparison between the FEniCSx focused ultrasound (FUS) solver and the \(k\)-Wave solver, where the source transducer is a focused bowl transducer. The simulation depicts a two-dimensional slice of a three-dimensional simulation on a subset of a human skull. This comparison is based on benchmark simulations from [1]. Figure courtesy of the author.
Figure 3. Pressure amplitude comparison between the FEniCSx focused ultrasound (FUS) solver and the \(k\)-Wave solver, where the source transducer is a focused bowl transducer. The simulation depicts a two-dimensional slice of a three-dimensional simulation on a subset of a human skull. This comparison is based on benchmark simulations from [1]. Figure courtesy of the author.

Our future work will incorporate nonlinearity into the mathematical model, which is vital to the study of high-intensity FUS applications. We also intend to explore the potential effects of graphics processing units and anticipate that graphics processors will be able to significantly improve solver performance, as vector assembly lends itself naturally to parallel computation.

Adeeb Arif Kor delivered a minisymposium presentation on this research at the 2023 SIAM Conference on Computational Science and Engineering (CSE23), which took place in Amsterdam, the Netherlands, last year. He received funding to attend CSE23 through a SIAM Student Travel Award. To learn more about Student Travel Awards and submit an application, visit the online page. 

SIAM Student Travel Awards are made possible in part by the generous support of our community. To make a gift to the Student Travel Fund, visit the SIAM website. 

References
[1] Aubry, J.-F., Bates, O., Boehm, C., Butts Pauly, K., Christensen, D., Cueto, C., … van ’t Wout, E. (2022). Benchmark problems for transcranial ultrasound simulation: Intercomparison of compressional wave models. J. Acoust. Soc. Am., 152, 1003-1019.
[2] Coussios, C.C., & Roy, R.A. (2008). Applications of acoustics and cavitation to noninvasive therapy and drug delivery. Annu. Rev. Fluid Mech., 40, 395-420.
[3] Duczek, S., & Gravenkamp, H. (2019). Mass lumping techniques in the spectral element method: On the equivalence of the row-sum, nodal quadrature, and diagonal scaling methods. Comput. Methods Appl. Mech. Eng., 353, 516-569.
[4] Lemprière, S. (2023). Ultrasound lets gene therapy into the brain. Nat. Rev. Neurol., 19(6), 326.
[5] Leung, S.A., Moore, D., Webb, T.D., Snell, J., Ghanouni, P., & Butts Pauly, K. (2021). Transcranial focused ultrasound phase correction using the hybrid angular spectrum method. Sci. Rep., 11(1), 6532.
[6] Meng, Y., Hynynen, K., & Lipsman, N. (2021). Applications of focused ultrasound in the brain: From thermoablation to drug delivery. Nat. Rev. Neurol., 17(1), 7-22.

About the Author

Adeeb Arif Kor

Ph.D. student, University of Cambridge

Adeeb Arif Kor is a Ph.D. student in the Department of Engineering at the University of Cambridge in the U.K. and a SLAB Fellow at the University of Malaya in Malaysia. His research interests include computational acoustics, high-performance computing, and open source software development.