Modeling ultrasound wave propagation

[bastpg]

Hi -- A colleague of mine, Mathias Davids, uses MFEM for electromagnetics modeling in the human body. We were impressed by the computation time and am wondering if this tool could be good for another biophysics problem of mine, i.e. modeling of propagation of ultrasound waves in the human body.

I am sure that MFEM can be used for this purpose, since many use the software for fluid dynamics, which is a more general and complex version of ultrasound wave propagation. The equations of the problem are nicely summarized here (Chapter 2, page 5): http://www.k-wave.org/manual/k-wave_user_manual_1.1.pdf

K-wave is a nice tool and could be used for validation but 1) it's slow and 2) it only supports voxel meshes. Does anyone have experience with this kind of physics problem in MFEM? I would love to hear thoughts on how to get started on this with minimal pain.

Thanks --

Bastien Guerin

[mlstowell]

Hello, @bastpg,

This looks like another fascinating biophysics application! Thank you for reaching out to us.

I glanced quickly at the manual you suggested. Implementing the coupled linear system in equations 2.1 should be relatively straightforward. Implementing the nonlinear equations described by equations 2.4 or 2.7 would, of course, be more difficult.

If equations 2.1 are sufficient for you to start experimenting with MFEM then I would suggest looking at the electromagnetics mini app Maxwell. This is also a coupled first order system but it models the electric and magnetic fields E and B. To modify this to model equations 2.1 I would suggest replacing E with the momentum rho_0 u and B with the density rho. The next modification would be to use the Raviart-Thomas (RT) basis functions in place of the Nedelec (ND), which are currently used for E, to represent the momentum. Also, change the Raviart-Thomas basis functions used for B to the discontinuous L2 basis functions for the density rho. Then you would need to change the differential operators replacing the curl with the divergence, and the weak curl with the weak gradient, etc..

I don't think Maxwell has a perfectly matched layer boundary condition but it shouldn't be too hard to add assuming your outer boundary isn't too complicated.

If the nonlinear equations are critical then I'd need to defer to other members of the team who are more familiar with non-linear PDEs (perhaps @jamiebramwell could offer her insights).

This may sound daunting but I'd be happy to walk you through the process whenever you run into trouble.

Best wishes,
Mark

[bastpg]

Hi Mark -- thanks for this thoughtful response. This sounds reasonable.

I am actually working with a slightly different formulation of the problem now, where the non-linearities are all approximated linearly. it's much simpler... This results in the following system of equations:
[Eq. 1]

Combining these equations yields a classic propagation equation with attenuation:
Eq. [2]

where the apparent attenuation coefficient is given by

and includes both attenuation and dispersion terms.

I have a couple of questions about this:

1. From the perspective of MFEM, is it better/easier to solve Eq. [1] or Eq. [2] (then followed by potentially and estimation of v using the second equation of the system Eq. [1])?
2. As an FEM technique, I expect MFEM to compute the steady state for the variable P and v, correct? That's fine, this is all I care about. I am not particularly interested in visualization of wave-front propagation or anything like that.
3. I know this will be almost impossible to answer. With that in mind... I am currently solving Eq. [1] using FDTD accelerated on a GPU. The FDTD method has a regular voxel grid, which is a downside. Another downside is that for an accurate computation, (a) the spatial resolution needs to be small (because of the finite difference scheme) and (b) the FDTD time-step needs to be small (because of the CFL stability criterion). Can you talk a bit about how these limitations compare to FEM? Is it reasonable to expect that the FEM problem is better conditioned than FDTD and thus allows a reduction of the spatial resolution?

Thanks,

Bastien Guerin

[mlstowell]

Hi, Bastien,

Regarding your first question, my personal preference would be to use equation 1. If you are interested in a time domain solution then the coupled first order equations are easier to integrate in time. MFEM has a symplectic integrator which should conserve energy nicely for this system of equations. You could also model the second order equation in the time domain I find the time integration to be less elegant and not necessarily energy preserving but it can be done in MFEM.

If, on the other hand, you are only interested in the steady state solution for an excitation at a specific frequency then the frequency domain would be a better choice. In this case equation 2 would be preferable, although again either equation could be used. The choice would be based on the available solvers and I suspect a block-diagonal multigrid preconditioner would work quite well with the second equation. I think that addresses question 2.

Question 3 is indeed more difficult. FEM methods do not avoid the CFL condition. The largest stable time step is computed from the largest eigenvalue of the time advance operator which, roughly speaking, is related to the smallest wavelength that can be represented by the basis functions. Rather than being determined by the mesh spacing, as in an FDTD method, this is determined by the spatial oscillations of the basis functions. If you use first order, i.e. piecewise linear, basis functions this leads to essentially the same CFL condition as FDTD. If you use piecewise quadratic or higher polynomial basis functions the CFL condition would be more restrictive than FDTD on the same mesh. However, higher order basis functions often allow you to use a coarser mesh which reduces the CFL condition. In either case the trick to simulating wave equations in either FDTD or FEM is to use a representation which accurately models the smallest wavelengths of interest and no smaller unless you're willing to pay the price of using a very small time step.

Off the top of my head I can think of three potential advantages to using an FEM method rather than finite difference and they are all related to geometric fidelity. If your problem domain is not brick-shaped or a relatively simple collection of bricks you may end up with many more degrees of freedom than are truly necessary to model the domain. Next, if your materials are inhomogeneous and hence produce waves of very different wavelengths in different regions then a finite difference grid will probably over resolve the longer wavelengths which again leads to more degrees of freedom than necessary. Finally, if your material interfaces are not aligned with the finite difference grid those interfaces may not be adequately captured. This, of course, assumes you can create a better aligned finite element mesh which may not be the case. I know such interface problems can lead to inaccuracies in eigenvalue computations but I can only guess similar problems may arise in computations of reflection/transmission coefficients.

As a final note I would add that finite difference methods might be better if the geometry involves many small features with sharp angles but relatively uniform wavelength. If the geometry contains larger, relatively homogeneous regions with occasional small sharp features I would expect FEM, with high order basis functions, to be the better choice.

Best wishes,
Mark