k-Wave User Forum » Topic: Scalloping in piston transducer simulation (large grid)

k-Wave User Forum » Topic: Scalloping in piston transducer simulation (large grid) http://www.k-wave.org/forum/topic/scalloping-in-piston-transducer-simulation-large-grid Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 18:27:33 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Scalloping in piston transducer simulation (large grid)" http://www.k-wave.org/forum/topic/scalloping-in-piston-transducer-simulation-large-grid#post-6438 Thu, 26 Apr 2018 09:13:45 +0000 Bradley Treeby 6438@http://www.k-wave.org/forum/ Hi apigula, The problem is related to the sampling of the time domain waveform. Although it might be sufficiently sampled (its Fourier spectrum decays), a time sample might not lie exactly on the maximum. A better way is to use exactly an integer number of points per period, and then extract the amplitude using the (e.g., <code>extractAmpPhase</code>). An alternative is to upsample the time domain waveforms (e.g., using <code>interpft</code>) before taking the maximum. Hope that helps, Brad. apigula on "Scalloping in piston transducer simulation (large grid)" http://www.k-wave.org/forum/topic/scalloping-in-piston-transducer-simulation-large-grid#post-6351 Mon, 19 Mar 2018 15:27:32 +0000 apigula 6351@http://www.k-wave.org/forum/ Yes, I increased it to allow a wave to propagate across the longest diagonal at speed C0, plus a little extra. Bradley Treeby on "Scalloping in piston transducer simulation (large grid)" http://www.k-wave.org/forum/topic/scalloping-in-piston-transducer-simulation-large-grid#post-6324 Wed, 07 Mar 2018 22:37:44 +0000 Bradley Treeby 6324@http://www.k-wave.org/forum/ I haven't run your code, but a quick question, did you also increase <code>T_END</code> when increasing the grid size? Brad. apigula on "Scalloping in piston transducer simulation (large grid)" http://www.k-wave.org/forum/topic/scalloping-in-piston-transducer-simulation-large-grid#post-6312 Mon, 05 Mar 2018 20:00:41 +0000 apigula 6312@http://www.k-wave.org/forum/ 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'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%