Validating the simulation results produced by Field II and FOCUS

When developing new simulation software, it is very helpful to have a known reference to check the results. When comparing new code to existing software, it is important to know whether the existing software is correct or not. We make an effort to validate our code whenever possible. We usually compare our results to the Raleigh-Sommerfeld integral, and we also evaluate comparisons with the impulse response and Field II. We have found that our results consistently converge to the same result as the Rayleigh-Sommerfeld integral, albeit much more rapidly. In fact, we have found that our method (namely the fast nearfield method or FNM) is a much better reference due to much faster convergence and much shorter computation times, but we still like to compare to something else just to make sure that the results are correct.

In these comparisons, we have uncovered some interesting behavior of Field II. In the text and figures below, we show what happens, and how to address these problems.

How to Model Attenuation of Time-Harmonic Pressures with Field II

This issue was originally described on Duo Chen’s web page: http://www.egr.msu.edu/~chenduo/. Duo found that Field II correctly attenuated the pressure amplitude for time-harmonic pressure calculations, but the phase obtained with Field II was incorrect. Duo provides Matlab scripts on his web page demonstrating this result. Duo showed that the frequency-dependent attenuation term is responsible for this incorrect phase. We have replicated these results in FOCUS, and the pressures are compared to those obtained with Field II.

The Field II manual states that attenuation should be implemented as follows:

set_field ('att', atten *100);
set_field ('freq_att', atten *100/f0);
set_field ('att_f0', f0);
set_field ('use_att', 1);

However, for time-harmonic pressure calculations, one of the values given above is incorrect. In the time-harmonic case, there is no need for frequency-dependent attenuation, and this term must be set to zero in order to obtain the correct phase. Thus, the correct settings are:

set_field ('att', atten*100);
set_field ('freq_att', 0);
set_field ('att_f0', f0);
set_field ('use_att', 1);

Results

Figure 1 shows the on-axis pressure results computed for a circular piston with radius = 4l excited by a time-harmonic signal with a frequency of 2MHz. Results are computed with the FNM and Rayleigh Sommerfeld integral routines provided by FOCUS and with the Field II program before and after the correction indicated above. The results show that the amplitudes are all in agreement.

Figure 2 shows the on-axis errors computed for the Rayleigh Sommerfeld integral, Field II before the attenuation term is corrected, and Field II after the attenuation term is corrected. The FNM expression provides the reference for these error calculations. Selecting the FNM as the reference is motivated by a comparison with the analytical expression for the on axis pressure (see Kinsler and Frey for this formula), which indicates that the FNM expression is accurate to 14 digits for the combination of parameters evaluated here. The error is computed by evaluating the difference between the complex pressure quantity in question and the complex reference pressure, taking the absolute value, and then dividing by the maximum absolute value of the reference pressure. This error metric includes the effects of phase and amplitude, or equivalently, errors in the real and imaginary parts of the pressure. As seen in Figure 2, the error in the complex pressure computed with Field II is substantial (with a peak value of about 36%) before the attenuation term is corrected, and the error is significantly smaller after the attenuation term is corrected.

Figure 3 shows the unwrapped phase error for Field II before and after the attenuation term is corrected and for the Rayleigh Sommerfeld integral. As before, the FNM provides the reference for the on axis pressure that is accurate to 14 digits. In Figure 3, the phase error is approximately equal to zero for the Rayleigh Sommerfeld integral and the Field II calculation with the corrected attenuation term. The result obtained with Field II before the attenuation term is corrected shows that the unwrapped phase error is linear with respect to distance when evaluated on axis for this set of parameters.

Implications

Such evaluations are important for us as we evaluate the FOCUS code and continue the development of new programs. For most applications, simulations that demonstrate such phase errors can be acceptable if we are only interested in the pressure amplitude or the approximate intensity. However, we are interested in knowing which answer is correct, and these comparisons (as well as comparisons with the analytical solution) provide us with that information. Furthermore, this could be an issue with time domain calculations, where phase errors are manifested as time shifts. We have not observed such errors in the transient results in Field II when attenuation is used, but we will look into this further at some point in the future.

Comparisons of the computation times and errors obtained with each method are also useful. We found that the FNM algorithm in FOCUS was more than 1000 times faster than Field II for this particular example. Furthermore, the FNM algorithm in FOCUS was accurate to 14 digits, whereas the corrected Field II program was accurate to only 3 digits. Links to scripts that repeat these results are provided so that individual users can make these comparisons themselves. To run these, make sure Field II and FOCUS are in the Matlab path, then execute the script compareFieldIItoFOCUS.m . Depending on the parameters of the problem, we typically find at least one or two orders of magnitude reduction in the computation time with FOCUS.

Other interesting behavior in Field II

Be careful turning attenuation on and off in Field II. The difference between

set_field ('use_att', 1);

and

set_field ('use_att', 0);

is significant, even if all of the other attenuation parameters are zero. Turning attenuation off multiplies all results by the sampling frequency f_s. If you scale your results after you compute them, you might not notice this. A better way to turn attenuation off is to always set ‘use_att’ to ‘1’ and then

set_field ('att', 0);

If you want to turn attenuation ‘on’ again, then insert the desired attenuation value in the line above.

Figure 1: On axis pressure amplitudes computed for a circular piston on-axis with the FNM and Rayleigh Sommerfeld integral routines provided by FOCUS and the Field II program before and after the correction of the attenuation factor. The computed amplitudes are all in agreement.

Figure 2: The normalized difference between the complex pressures evaluated on axis where the FNM provides a reference that is accurate to 14 digits. This normalized error is computed for Field II before the attenuation factor is corrected, Field II after the attenuation term is corrected, and the Rayleigh Sommerfeld integral. The errors are large before the attenuation term is corrected (reaching a maximum of about 36%), and the errors are much smaller after the attenuation term is corrected and when the pressure is calculated with the Rayleigh Sommerfeld integral.

Figure 3: Unwrapped phase errors calculated with Field II before the attenuation term is corrected, with Field II after the attenuation term is corrected, and with the Rayleigh Sommerfeld integral. The phase errors are close to zero with Field II after the attenuation term is corrected and with the Rayleigh Sommerfeld integral. The phase error increases linearly with distance with the erroneous attenuation term in Field II.

MATLAB Attachments

Contributors

Simulations

Web Site