Hi,
I'm currently trying to simulate the attenuation in a homogeneous medium (it happens to be bone) in 3D. My source is a simple disc lying centred in the one of the grid planes. I'm trying to confirm that I do indeed get the expected attenuation that I might expect from the values that apply to the medium. The properties are:
medium.alpha_coeff = 13.02;
medium.alpha_power = 1.03;
medium.speed_of_sound = 3000;
medium.density = 1850;
The power is set to 1.03, as I want it to essentially be inversely linear with frequency, as a simple approximation.
The simulation is being conducted with a CW of 0.5MHz. There is a PML of 10 pts and an alpha of 2 at each boundary. The grid spacing is 0.5mm in all dimensions, and the stepsize is 0.5e-8s.
I take a measurement along the axial direction of the beam created by the disc, and either take RMS values of the last X cycles of the signal at each location, or the maximum value reached throughout the time series at each location. It's my understanding that once out of the near field, the wave along this axis will essentially be planar.
I have already run 1D simulations of a similar scenario, and once turning dispersion off, saw the behaviour that I would expect. However, in 3D I just can't seem to get the attenuation as predicted (~6.5dB/cm @ 0.5MHz). In addition, I seem to see some odd effects. If the alpha_coeff is any value greater than zero, the wavelength of the CW seems to increase greater than I would expect. This I certainly did not expect and has confused me greatly. Interestingly, if I use the source.p0 approach then the delta function seems to decay as expected. In addition, when running the 3D model turning dispersion on or off makes no difference, unlike in the 1D case.
Below I've include the 3D simulation code.
Thanks in advance for any help,
Regards,
Lewis
% --------------------------------------------------------------------
clear;
%% Define grid.
Nx = 128;
Ny = 512;
Nz = 128;
dx = 0.5e-3;
dy = 0.5e-3;
dz = 0.5e-3;
kgrid = kWaveGrid(Nx,dx,Ny,dy,Nz,dz);
kgrid.dt = 0.5e-8;
kgrid.Nt = 5000;
pml_size = [10 10 10];
pml_alpha = 2;
%% Define medium.
medium.alpha_coeff = ones(Nx,Ny,Nz)*13.02;
medium.alpha_power = 1.03; % Must be a scalar.
medium.sound_speed = ones(Nx,Ny,Nz)*3000;
medium.density = ones(Nx,Ny,Nz)*1850;
%% Define source.
src_type = 'circ_piston';
source.p_mask = false(Nx,Ny,Nz);
switch src_type
case 'monopole'
source.p_mask(Nx/2,Ny/2,Nz/2) = true;
case 'circ_piston'
disc = makeDisc(Nx,Nz,0,0,Nx/4);
source.p_mask(:,Ny/2,:) = disc;
end
src_freq_hz = 500000;
source.p = sin(2*pi*src_freq_hz*kgrid.t_array);
%% Define Sensor positions.
sensor.mask = false(Nx,Ny,Nz);
sensor.mask(:,:,Nz/2) = true;
%% Run CUDA simulation.
input_args = {...
'PlotLayout',false,'PlotPML',false,...
'PMLInside',true,'PMLSize',pml_size,'PMLAlpha',pml_alpha...
'DataCast', 'single', 'CartInterp', 'nearest'};
outputPath = fullfile(def_default_dir('kwave_data'),'kw_test'); % Will need to be changed.
input_args = [input_args ...
{'SaveToDisk',outputPath}];
kspaceFirstOrder3D(kgrid,medium,source,sensor,input_args{:});
cuda_input_args = {...
'BinaryPath' fullfile(def_default_dir('kwave_bin'),'cuda'),... % Will need to be changed.
'DataPath' def_default_dir('kwave_data'),... % Will need to be changed.
'DataName' 'kw_test',...
'DeleteData' false,...
'PMLInside' true,...
'PMLSize' pml_size,...
'PMLAlpha' pml_alpha,...
'DataCast' 'single',...
'CartInterp' 'nearest',...
'BinaryName' 'kspaceFirstOrder3D-CUDA.exe',...
'DeviceNum' 0};
sensor_data = kspaceFirstOrder3DG(kgrid,medium,source,sensor,cuda_input_args{:});
%% Plot results.
% Get sensor pressure data from sim.
output_filename = fullfile(def_default_dir('kwave_data'),'kw_test_output.h5');
sensor_output = h5read(output_filename,'/p',[1 1 1],[inf inf inf]);
p = reshape(sensor_output,[Nx Ny kgrid.Nt]);
p_axial = squeeze(p(Nx/2,:,:));
figure;
plot((1:Ny)*dy*1000,20*log10(max(p_axial,[],2)));
ylabel('spl (dB)')
xlabel('axial dist. (mm)')
figure;
plot((1:Ny)*dy*1000,p_axial(:,end));
ylabel('amplitude')
xlabel('axial dist. (mm)')
%% END.