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');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%