Hello, I have two questions about large-grid simulations, based on a modified version of your circular piston example.
(1) At large grid sizes (say, 256x256x3456), I was getting terrible agreement between the theory and simulated axial pressures in the farfield. I was able to improve the simulation's accuracy significantly by using PML = 30. This is nonintuitive to me - Why would a large grid size require a larger PML?
(2) Even after increasing the PML, I get funny oscillations in the simulated axial pressure when using a large grid size. In particular, there's a scallop-like pattern near the farfield peak, and the pressure still doesn't quite match the prediction. Any idea why this is?
I'd appreciate any advice!
Here's the code I was running:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
% =========================================================================
% DEFINE LITERALS
% =========================================================================
% select which code to run
% 1: MATLAB CPU code
% 2: MATLAB GPU code
% 3: C++ code
MODEL = 3;
% scale parameter (changes the resolution of the simulation)
SC = 1;
% grid parameters
PML_SIZE = 30; % size of the perfectly matched layer [grid points]
Nx = 3456*SC - 2*PML_SIZE; % number of grid points in the x direction
Ny = 256*SC - 2*PML_SIZE; % number of grid points in the y direction
Nz = 256*SC - 2*PML_SIZE; % number of grid points in the z direction
% sampling parameters
PPW = 5*SC; % points per wavelength (this defines the grid size)
PPP = 80*SC/2; % points per period (this defines the CFL)
T_END = 1.3e-4; % total compute time [s] (this must be long enough to reach steady state)
RECORD_PERIODS = 2; % number of periods to record
% medium parameters
C0 = 1580; % sound speed [m/s]
RHO0 = 1000; % density [kg/m^3]
% source parameters
F0 = 6e6; % source frequency [Hz]
SOURCE_RADIUS = 90*SC; % piston radius [grid points]
SOURCE_MAG = 0.5e6; % source pressure [Pa]
% =========================================================================
% CREATE GRID
% =========================================================================
% calculate the grid spacing based on the PPW and F0
dx = C0 / (PPW * F0); % [m]
dt = 1 / (PPP*F0); % [s]
% calculate the CFL
disp(['CFL = ' num2str(C0*dt/dx)]);
% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx);
% create the time array
kgrid.t_array = 0:dt:T_END;
% assign medium properties
medium.sound_speed = C0;
medium.density = RHO0;
% =========================================================================
% CREATE SOURCE
% =========================================================================
% create piston source mask
piston = makeDisc(Ny, Nz, Ny/2, Nz/2, SOURCE_RADIUS);
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(1, :, :) = piston;
% create time varying source
source.p = SOURCE_MAG * sin(2*pi*F0*kgrid.t_array);
% filter source with up-ramp
ramp_length = round((2*pi/F0)/dt);
ramp = (-cos( (0:(ramp_length-1))*pi/ramp_length ) + 1)/2;
source.p(1:ramp_length) = ramp .* source.p(1:ramp_length);
% =========================================================================
% CREATE SENSOR MASK
% =========================================================================
% set sensor mask to record central axis, not including the source point
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(2:end, Ny/2, Nz/2) = 1;
% record the maximum pressure
sensor.record = {'p_max'};
% average only the final few periods when the field is in steady state
sensor.record_start_index = kgrid.Nt - RECORD_PERIODS*PPP + 1;
% =========================================================================
% RUN k-WAVE SIMULATION
% =========================================================================
% run code
switch MODEL
case 1
% MATLAB CPU
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...
'DataCast', 'single', 'PMLInside', false, ...
'DisplayMask', source.p_mask, 'PlotScale', [-1, 1]*SOURCE_MAG);
case 2
% MATLAB GPU
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...
'DataCast', 'gpuArray-single', 'PMLInside', false, ...
'DisplayMask', source.p_mask, 'PlotScale', [-1, 1]*SOURCE_MAG);
case 3
% C++
sensor_data = kspaceFirstOrder3DC(kgrid, medium, source, sensor, 'PMLInside', false, 'PMLSize', PML_SIZE);
end
% =========================================================================
% ANALYTICAL SOLUTION
% =========================================================================
% calculate the wavenumber
k = 2*pi*F0./C0;
% define radius and axis
a = SOURCE_RADIUS*dx;
x = (kgrid.x_vec(2:end, :) - min(kgrid.x_vec(:)));
% calculate the analytical solution for a piston in an infinite baffle
% for comparison (Eq 5-7.3 in Pierce)
r_a = sqrt(x.^2 + a^2);
p_ref = SOURCE_MAG*abs(-2i*exp(1i*k*(x+r_a)/2).*sin((k*r_a - k*x)/2));
% =========================================================================
% VISUALISATION
% =========================================================================
% plot the pressure along the focal axis of the piston
figure;
plot(1000*x, p_ref, 'k-');
hold on;
plot(1000*x, sensor_data.p_max, 'bx')
xlabel('Distance (cm)');
legend('Exact', 'k-Wave');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%