Dear Treeby and Colleagues,
I am attempting to run a convergence test for a simple simulation, but am not seeing convergence. I am simulating a single 9 mm element at 650 kHz, and am concerned with the pressure field on axis 15 cm away (this is to simulate the element's contribution in a larger array). My results are this:
Voxel size (mm) || Pressure(@15 cm)
0.2 .. 365
0.3 .. 78.9
0.4 .. 266
0.5 .. 199
For the convergence test, I am keeping the grid size constant, keeping the actual size of the transducer (a flat single element) constant by changing its number of pixels, and updating the timestep and continuous wave source using makeTime() and createCWSignals(). The PML does not reach the regions of interest in any of the grid configurations. One thing that stands out as unique is that my grid is [64 64 850] spatial points so is rather long and skinny.
Does anything stand out to you as notable as to why the pressure in the far field is changing dramatically?
Thank you!
Here is my code:
% script to simulate a single flat element
clear all
f0 = 0.65E6; % frequency [MHz]
Amp = 1000; % amplitude at source.p_mask
phase = 0; % phase of source wave
sensor.record = {'p_max'}; %set to p_max, or p_rms depending on what you want
%% pick initial dimensions
vox = .3E-3;
dim = [60 60 850];
'Your simpace size'
dim
%checkFactors(min(dim),max(dim)+100)
%% look at the dimensions of your image and select a final grid size that fits
% your simspace and is in one of the first 3 rows output by checkFactors
padsize = [64 64 850];
%% setup and run simulation
% pad arrays
foo = ones(padsize);
dim = padsize;
% get medium properties (density and speed of sound for now)
medium.sound_speed = 1480.*foo;
medium.density = 997.*foo;
% create kgrid
kgrid = kWaveGrid(dim(1), vox, dim(2), vox, dim(3), vox);
[kgrid.t_array, dt_SIM] = makeTime(kgrid, medium.sound_speed);
% create pressure vector
source.p = createCWSignals(kgrid.t_array, f0, Amp, phase);
%% make transducer
foc = [32 32 800]; % location at region of interest (focus of larger array)
roc = 150; %[mm] % radius of curvature of larger array
D = round((9E-3)/vox)+1; %[m] element diameter
rocp = round((roc.*1E-3)/vox);
top = [foc(1) foc(2) foc(3)-rocp]; % place the element 150 mm away from foc
%I want to simulate a flat element, so I change the radius to be huge
rocp = 10000;
xdcr = makeBowl(dim,top,rocp,D,foc);
%%
source.p_mask = xdcr==1; %only get spherical cap part of xdcr file
sensor.mask = ones(dim);
%%
sensor_data = kspaceFirstOrder3DG(kgrid,medium,source,sensor);
% grab data
if sensor.record{1} == 'p_rms'
data = sensor_data.p_rms;
elseif sensor.record{1} == 'p_max'
data = sensor_data.p_max;
else
'this script only supports recording p and p_rms'
end
%
data = reshape(data,dim);