Hi,
I modified modeling non-linearity example to run it on the air at 40 kHz. I have different results for k-Wave and Mendousse simulations. I set the wave separation to be less than shock formation distance. The file looks like this
*****
%
% define the properties used in the simulation
p0 = 290; % source pressure [Pa]
c0 = 343; % sound speed [m/s]
rho0 = 1.225; % density [kg/m^3]
alpha_0 = 0.25; % absorption coefficient [dB/(MHz^2 cm)]
% sigma = 2; % shock parameter
source_freq = 40000; % frequency [Hz]
points_per_wavelength = 100; % number of grid points per wavelength at f0
% wavelength_separation = 10; % separation between the source and detector
wavelength_separation = 40;
pml_size = 80; % PML size
pml_alpha = 1.5; % PML absorption coefficient [Np/grid point]
CFL = 0.25; % CFL number
BonA = .7;
% =========================================================================
% RUN SIMULATION
% =========================================================================
% compute grid spacing
dx = c0 / (points_per_wavelength * source_freq); % [m]
% compute grid size
Nx = wavelength_separation * points_per_wavelength + 20; % [grid points]
% create the computational grid
kgrid = kWaveGrid(Nx, dx);
% assign the properties of the propagation medium
medium.sound_speed = c0;
medium.density = rho0;
medium.alpha_power = 2;
medium.alpha_coeff = alpha_0;
% extract the maximum frequency supported by the grid
f_max = min(medium.sound_speed) / (2 * dx); % [Hz]
% define a single source element
source_pos = 1;
source.p_mask = zeros(Nx, 1);
source.p_mask(source_pos) = 1;
% define a single sensor position an integer number of wavelengths away
sensor.mask = zeros(Nx, 1);
x_px = wavelength_separation * points_per_wavelength;
sensor.mask(source_pos + x_px) = 1;
x = x_px * dx
% shock formation distance rho0*c0^3/(omega*beta*p0)
xsh = rho0*c0^3/(2*pi*source_freq*(1+BonA)*p0)
medium.BonA = BonA;
% set the simulation options
input_args = {'PlotFreq', 80, 'PlotScale', [-1, 1] * p0 * 1.05, ...
'PMLInside', false, 'PMLSize', pml_size, 'PMLAlpha', pml_alpha};
% compute points per temporal period
points_per_period = round(points_per_wavelength / CFL);
% compute corresponding time spacing
dt = 1 / (points_per_period * source_freq);
% create the time array using an integer number of points per period
% t_end = 25e-6; %original
t_end = 2e-3;
Nt = round(t_end / dt);
kgrid.setTime(Nt, dt);
% set the start time to only record the last three periods
sensor.record_start_index = kgrid.Nt - 5 * points_per_period + 1;
% create the source term
source.p = p0 * sin(2 * pi * source_freq * kgrid.t_array);
% run the simulation
sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, input_args{:});
% create time axis for mendousse solution
t_axis = (0:(length(sensor_data) - 1)) .* dt;
% compute mendousse solution for comparison
p_mendousse = mendousse(x * ones(size(t_axis)), t_axis, source_freq, p0, c0, rho0, BonA, alpha_0);
% get the amplitude spectra
as_mendousse = spect(p_mendousse, 1/dt);
% f = f * 1e-6;
as_kspace = spect(sensor_data, 1/dt);
% =========================================================================
% VISUALISATION
% =========================================================================
% plot the time series
figure;
subplot(3, 1, 1);
plot(t_axis, sensor_data, 'r-');
hold on;
plot(t_axis, p_mendousse, 'k--');
xlabel('Time');
ylabel('Pressure');
legend('k-Wave', 'Mendousse');
subplot(3, 1, 2);
% plot simulation
plot(as_kspace,'color','r');
legend('k-Wave');
subplot(3, 1, 3);
% plot reference
plot(as_mendousse, 'color', 'k');
legend('Mendousse');
*******
thank you very much!
vadim