k-Wave User Forum » Topic: Simulation Unstable

k-Wave User Forum » Topic: Simulation Unstable http://www.k-wave.org/forum/topic/simulation-unstable Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 21:36:17 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php ArteryLee on "Simulation Unstable" http://www.k-wave.org/forum/topic/simulation-unstable#post-9128 Wed, 21 Aug 2024 16:47:26 +0000 ArteryLee 9128@http://www.k-wave.org/forum/ Hi Bradley, Thanks very much for your reply. I have tried to remove the line <code>medium.alpha_mode = 'no_dispersion'</code>. However, it still produces all NaN matrix in the sensor in this case as long as I assign alpha coefficient to the collimator mask by <code>medium.alpha_coeff(collimator_mask) = alpha_collimator</code>. I mean, once I comment out this line, it can produce non-NaN results with the C++/CUDA code. BTW, I just modified my code according to your advice by setting the power law exponent to 2. In this case, it can also produce non-NaN results with the C++/CUDA code. But I am not sure how to get the correct power law absorption coefficients based on my previous numbers now. Would you mind explaining more about this? Thank you! Best, Artery Bradley Treeby on "Simulation Unstable" http://www.k-wave.org/forum/topic/simulation-unstable#post-9126 Sun, 18 Aug 2024 10:36:50 +0000 Bradley Treeby 9126@http://www.k-wave.org/forum/ Hi Artery, the reason is that the no dispersion input is not supported in the C++ code. Unfortunately the docs are not super clear on this. However, if you’re doing a CW simulation, you can also get no dispersion but setting the power law exponent to 2, and then adjusting the pre factor so you get the same level of absorption at your driving frequency. ArteryLee on "Simulation Unstable" http://www.k-wave.org/forum/topic/simulation-unstable#post-9123 Wed, 14 Aug 2024 19:39:54 +0000 ArteryLee 9123@http://www.k-wave.org/forum/ Hi Bradley, I encountered a situation similar to the one Toni described above. In my simulation, I added a collimator into the medium. With <code>kspaceFirstOrder3DC</code> or <code>kspaceFirstOrder3DG</code>, I got all NaN in the sensor grid as long as I added the collimator. While with common <code>kspaceFirstOrder3D</code>, whether casting data to <code>single</code> or <code>gpuArray-single</code>, the results were normal. Meanwhile, without adding the collimator, <code>kspaceFirstOrder3DC</code> and <code>kspaceFirstOrder3DG</code> work fine with pure water medium. I'm not sure how this would happen as I have used the function <code>checkStability</code> to make sure my dt (3.3333e-08) is smaller than the dt_stability_limit (5.0456e-08). Thank you very much for taking a look into this! Best, Artery %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ========================================================================= % DEFINE LITERALS % ========================================================================= % medium parameters c0 = 1500; % sound speed [m/s] rho0 = 1000; % density [kg/m^3] alpha0 = 3.0227e-05; % water power law absorption coefficient [dB/(MHz^y cm)] % collimator parameters c_collimator = 2410; rho_collimator = 1160; alpha_collimator= 9.54; % source parameters source_f0 = 5e5; % source frequency [Hz] source_roc = 59.3e-3; % bowl radius of curvature [m] source_diameter = 25.4e-3; % bowl aperture diameter [m] source_amp = 1.4*2e4; % source pressure [Pa] focus_depth = 38.1e-3; % [m] % grid parameters axial_size = 70e-3; % total grid size in the axial dimension [m] lateral_size = 70e-3; % total grid size in the lateral dimension [m] % computational parameters ppw = 12; % number of points per wavelength t_end = 200e-6; % total compute time [s] record_periods = 1; % number of periods to record cfl = 0.2; % CFL number source_x_offset = 20; % grid points to offset the source bli_tolerance = 0.01; % tolerance for truncation of the off-grid source points upsampling_rate = 10; % density of integration points relative to grid % ========================================================================= % DEFINE KGRID % ========================================================================= % calculate the grid spacing based on the PPW and F0 dx = c0 / (ppw * source_f0); % [m] % compute the size of the grid Nx = roundEven(axial_size / dx) + source_x_offset; Ny = roundEven(lateral_size / dx); Nz = Ny; % create the computational grid kgrid = kWaveGrid(Nx, dx, Ny, dx, Nz, dx); % compute points per temporal period PPP = round(ppw / cfl); % compute corresponding time spacing dt = 1 / (PPP * source_f0); % create the time array using an integer number of points per period Nt = round(t_end / dt); kgrid.setTime(Nt, dt); % ========================================================================= % DEFINE KWAVEARRAY FOR TRANDUCER % ========================================================================= % set bowl position and orientation bowl_pos = [kgrid.x_vec(1) + source_x_offset * kgrid.dx, 0, 0]; focus_pos = [bowl_pos(1) + focus_depth, 0, 0]; % create empty kWaveArray karray = kWaveArray('BLITolerance', bli_tolerance, 'UpsamplingRate', upsampling_rate, 'SinglePrecision', true); % add bowl shaped element karray.addBowlElement(bowl_pos, source_roc, source_diameter, focus_pos); % ========================================================================= % DEFINE SOURCE SIGNAL % ========================================================================= % create time varying source source_sig = createCWSignals(kgrid.t_array, source_f0, source_amp, 0); % assign binary mask source.p_mask = karray.getArrayBinaryMask(kgrid); % assign source signals source.p = karray.getDistributedSourceSignal(kgrid, source_sig); % ========================================================================= % DEFINE BINARY COLLIMATOR MASK % ========================================================================= [OUTPUTgrid] = VOXELISE(166,166,152,'Collimator_14mm.STL','xyz'); OUTPUTgrid = imrotate3(double(OUTPUTgrid), 90, [1 0 0]); OUTPUTgrid = logical(OUTPUTgrid); collimator_mask = zeros(Nx, Ny, Nz, 'logical'); collimator_mask(ceil(Nx/2)-ceil(size(OUTPUTgrid,1)/2)-59-14:ceil(Nx/2)+ceil(size(OUTPUTgrid,1)/2)-1-59-14, ... ceil(Ny/2)-ceil(size(OUTPUTgrid,2)/2)+1:ceil(Ny/2)+ceil(size(OUTPUTgrid,2)/2)-1+1, ... ceil(Nz/2)-ceil(size(OUTPUTgrid,3)/2)+2:ceil(Nz/2)+ceil(size(OUTPUTgrid,3)/2)-1+2) = OUTPUTgrid; % ========================================================================= % DEFINE MEDIUM % ========================================================================= % assign medium properties medium.sound_speed = c0.*ones(Nx, Ny, Nz); medium.sound_speed(collimator_mask) = c_collimator; medium.density = rho0.*ones(Nx, Ny, Nz); medium.density(collimator_mask) = rho_collimator; medium.alpha_coeff = alpha0.*ones(Nx, Ny, Nz); medium.alpha_coeff(collimator_mask) = alpha_collimator; medium.alpha_power = 1.001; medium.alpha_mode = 'no_dispersion'; % ========================================================================= % DEFINE SENSOR % ========================================================================= % set sensor mask to record central plane, not including the source point sensor.mask = zeros(Nx, Ny, Nz); sensor.mask(source_x_offset + 2:end, :, Nz/2 + 1) = 1; % mid-plane/central plane % record the pressure sensor.record = {'p', 'p_rms'}; % record only the final few periods when the field is in steady state sensor.record_start_index = kgrid.Nt - record_periods * PPP + 1; % ========================================================================= % SIMULATION % ========================================================================= % set input options input_args = {... 'PMLSize', 'auto', ... 'PMLInside', false, ... 'PlotPML', false, ... 'DisplayMask', 'off'}; % run code sensor_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor, input_args{:}); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Bradley Treeby on "Simulation Unstable" http://www.k-wave.org/forum/topic/simulation-unstable#post-6908 Tue, 25 Jun 2019 06:04:38 +0000 Bradley Treeby 6908@http://www.k-wave.org/forum/ Hi Toni, You can use the function <code>checkStability</code> to get an estimate of the maximum stable time step, and then use this with the <code>setTime</code> method of the <code>kWaveGrid</code> class. Brad. guns2roses on "Simulation Unstable" http://www.k-wave.org/forum/topic/simulation-unstable#post-6886 Mon, 17 Jun 2019 14:30:25 +0000 guns2roses 6886@http://www.k-wave.org/forum/ Hi Brad, I recently ran a simulation using "kspaceFirstOrder3DG" with multiple media as shown below. However, the simulation seems unstable as the result of "sensor_data" is a matrix full of "NaN (Not-a-Number)" in Matlab. I understand the differences of the "speed of sound" and "absorption" among different media are quite large which may be the reason causing simulation instability. I wonder if there is a way to tweak the values of "speed of sound" or "absorption" so that the simulation becomes stable while still maintaining certain amount of accuracy. Thanks a million! Best wishes, Toni %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% np2db=0.2*log10(exp(1)); ppts=zeros(5,5); % Velocity, Density, Absorption Coefficient, Power Law Exponent, Non-linearity (m/s, kg/m3, Np/cm/MHz, y, B/A) ppts(1,:)=[1561 1081 3.917*np2db 1.2598 6.756]; % heart muscle ppts(2,:)=[1482 994 0.025*np2db 1 4.96]; % water ppts(3,:)=[1578 1050 2.3676*np2db 1.05 6.1125]; % blood ppts(4,:)=[3515 1908 54.553*np2db 1 1]; % sternum ppts(5,:)=[1500 1050 161.88*np2db 1 1]; % lung % Please note that the 4th column "power law exponent" of matrix "ppts" has been disabled. % Instead the "medium.alpha_power" is set to 1.05 for all media as recommended in the k-Wave manual.