Hello,
I am running a series of simulations in which I have an array of pressure sources (all driven at the same amplitude and phase), with sensor points along the axis of the array. The entire simulation medium is air, with attenuation and nonlinearity included. The amplitude/envelope of the signal that is input to each pressure source is constant over time.
For a "small" grid size (320x320x128), the pressure at the sensor points is more or less what I expect. However, when I increase the grid size in the axial direction (320x320x1024), I find that the amplitude/envelope of the pressure at the sensor points increases with time, even though I've changed nothing else about the simulation.
I'm not sure what the reason for this might be - is this a symptom of some numerical instability? I don't see how this can be physical, but if it is, why does it only show up when I increase the grid size in the axial direction?
I've included the code below for reference. Thanks in advance for your help.
%%%%%%%% Setup %%%%%%%%
% create the computational grid
Nx = 1024; % number of grid points in the x direction
Ny = 320; % number of grid points in the y direction
Nz = 320; % number of grid points in the z direction
dx = 0.85e-3; % grid point spacing in the x direction [m]
dy = 0.85e-3; % grid point spacing in the y direction [m]
dz = 0.85e-3; % grid point spacing in the z direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
% define the properties of the propagation medium
clear('medium');
medium.sound_speed = 343; % [m/s]
medium.density = 1.225; % [kg/m^3]
medium.alpha_coeff = 0.8; % [dB/(MHz^y cm)]
medium.alpha_power = 1.33; % [y in fmla above]
medium.BonA = 0.4; % equal to (gamma - 1) for gases
% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
% define a square source element
source_width = 60; % [# of grid points]
source.p_mask = zeros(Nx, Ny, Nz);
source_x = 40;
source_y1 = 2*round(0.5*(Ny/2 - source_width/2 + 1));
source_y2 = 2*round(0.5*(Ny/2 + source_width/2));
source_z1 = 2*round(0.5*(Nz/2 - source_width/2 + 1));
source_z2 = 2*round(0.5*(Nz/2 + source_width/2));
for y = source_y1:4:source_y2
for z = source_z1:4:source_z2
source.p_mask(source_x,y,z) = 1;
end
end
% define the sensor mask
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(:, Ny/2, Nz/2) = 1;
% define a time-varying sinusoidal source
source_freq = 51.5e3; % [Hz]
source_mag = 3168; % [Pa]
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);
% smooth the source to avoid aliasing
source.p = filterTimeSeries(kgrid, medium, source.p,'ZeroPhase',true);
% enforce input pressure rather than adding it to the existing pressure
source.p_mode = 'additive';
% define the field parameters to record
sensor.record = {'p', 'u'};
% input arguments
input_args = {'DisplayMask', source.p_mask, 'DataCast', 'single', ...
'CartInterp', 'nearest', 'Smooth', false, 'PMLSize', 10, ...
'PMLAlpha', 2};
fprintf('\nSimulation parameters created.\n');
% view the source element and sensor mask
% voxelPlot( single( source.p_mask | sensor.mask ) );
%%%%%%%% Run %%%%%%%%
% Run directly in MATLAB
% sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
% Save input files to disk to be run with compiled C++ code
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...
input_args{:}, 'SaveToDisk', '5x5_nonlinear_prop_input.h5');
fprintf('\nData saved to HDF file.\n');
%%%%%%%% Load results %%%%%%%%
sensor_data_C.p = h5read('5x5_nonlinear_prop_output.h5', '/p');
sensor_data_C.ux = h5read('5x5_nonlinear_prop_output.h5', '/ux');
sensor_data_C.uy = h5read('5x5_nonlinear_prop_output.h5', '/uy');
sensor_data_C.uz = h5read('5x5_nonlinear_prop_output.h5', '/uz');
fprintf('\nData read from HDF file.\n');
fprintf('\nTime vector has %d entries.\n',size(sensor_data_C.p,2));
x = 0:(Nx-1);
x = x - source_x;
x = dx*x;
%%%%%%%% Visualize transient results %%%%%%%%
% pressure vs time at specific positions
figure
plot(sensor_data_C.p(60,:),'k');