k-Wave User Forum » Topic: difference in Mendousse and k-Wave simulation

k-Wave User Forum » Topic: difference in Mendousse and k-Wave simulation http://www.k-wave.org/forum/topic/difference-in-mendousse-and-k-wave-simulation Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:15:41 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php vpiskun on "difference in Mendousse and k-Wave simulation" http://www.k-wave.org/forum/topic/difference-in-mendousse-and-k-wave-simulation#post-6433 Fri, 20 Apr 2018 20:53:22 +0000 vpiskun 6433@http://www.k-wave.org/forum/ Thank you Brad! vadim Bradley Treeby on "difference in Mendousse and k-Wave simulation" http://www.k-wave.org/forum/topic/difference-in-mendousse-and-k-wave-simulation#post-6418 Wed, 18 Apr 2018 09:16:28 +0000 Bradley Treeby 6418@http://www.k-wave.org/forum/ I've never noticed this before, but it looks the <code>mendousse</code> function breaks for large values of the Gol'dberg number. If you reduce this (e.g., by increasing absorption to say <code>alpha0 = 4</code>), it works ok. If the shock parameter is also high, you could try coding up the Fay solution (described on page 136 of Nonlinear Acoustics, by Hamilton and Blackstock). Hope that helps, Brad. vpiskun on "difference in Mendousse and k-Wave simulation" http://www.k-wave.org/forum/topic/difference-in-mendousse-and-k-wave-simulation#post-6417 Wed, 18 Apr 2018 00:54:56 +0000 vpiskun 6417@http://www.k-wave.org/forum/ 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