k-Wave User Forum » Topic: Unexpected instability

k-Wave User Forum » Topic: Unexpected instability http://www.k-wave.org/forum/topic/unexpected-instability Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 06:07:24 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Foufour on "Unexpected instability" http://www.k-wave.org/forum/topic/unexpected-instability#post-6538 Fri, 27 Jul 2018 20:16:42 +0000 Foufour 6538@http://www.k-wave.org/forum/ Hi Brad, I didn't find a solution in the meantime. I even switched off smoothing, and smoothed manually. Same result. I replaced the external source file by a simple ball, and got the same instability again. So yes, I can create an example that doesn't depend on external functions or data that has the same problem. The only thing you'll have to change is the path of the created HDF5 input file (if you want to create one at all). The code follows... <pre><code>%% Input parameters grid_length = 0.1; % length of kgrid in m grid_width = 0.1; % width of kgrid in m axial_spacing = 250e-6; % spacing between neighbouring points in axial direction (dx) in m lateral_spacing = 250e-6; % spacing between neighbouring points in lateral direction (dy,dz) in m time_step = 25e-9; % simulation time step in s total_time = 35e-6; % total simulation time in s T = 25; % temperature in °C d_geant_grid = 100e-6; % grid spacing of the above source file smooth_p0 = true; % Boolean controlling whether initial pressure distribution gets smoothed automatically smooth_density = true; % Boolean controlling whether density distribution gets smoothed automatically smooth_soundspeed = true; % Boolean controlling whether sound speed distribution gets smoothed automatically StreamToDisk = true; % Boolean controlling whether simulation data are stored on disk during simulation DataCast = 'gpuArray-double'; % Variable precision & simulation mode. Possible: 'single', 'double', 'gpuArray-double', 'gpuArray-single' SaveToDisk = 'sirmiobeam-double-250-in'; % File name of the HDF5 input file PMLInside = false; %% Creating kWaveGrid % In all dimensions, the grid should have as small prime factors as % possible. That's why I chose to increase the grid to the next power of 2. % The 20 points subtracted will get re-added in the end as PML. Nx = 2^nextpow2(grid_length/axial_spacing)-20; Ny = 2^nextpow2(grid_width/lateral_spacing)-20; Nz = Ny; kgrid = kWaveGrid(Nx,axial_spacing,Ny,lateral_spacing,Nz,lateral_spacing); kgrid.dt = time_step; kgrid.Nt = round(total_time/time_step); %% Creating medium sound_speedwater = 1.40238677e3 + 5.03798765*T - 5.80980033e-2*T.^2 ... + 3.34296650e-4*T.^3 - 1.47936902e-6*T.^4 + 3.14893508e-9*T.^5; %[m/s] densitywater = 999.9 + T.*(0.06079 + T.*(-0.008313 + T.*(6.595e-5 + T.* ... (-4.437e-7 + 1.369e-9.*T))));%[kg/m^3] medium = struct('sound_speed', sound_speedwater,... 'density', densitywater, 'sound_speed_ref', sound_speedwater); %% Creating source raw_source.p0 = makeBall(500,1000,1000,200,500,500,50); % This is similar to what my function new_read_Geant does - it creates a % 500x1000x1000 source matrix. I simplified this to a ball source instead of % a Geant4 simulation output for you. The behaviour in terms of instablilty % is the same. %% Adjusting the source to the spacing of the k-Wave grid % This part corresponds to my adjust_grids function, in an attempt to % become independent of the inherent grid spacing of the source file. [x_rawlength, y_rawlength, z_rawlength] = size(raw_source.p0); x_newlength = ceil(x_rawlength * d_geant_grid / kgrid.dx); y_newlength = ceil(y_rawlength * d_geant_grid / kgrid.dy); z_newlength = ceil(z_rawlength * d_geant_grid / kgrid.dz); [y_geant, x_geant, z_geant] = meshgrid((1 : y_rawlength) * d_geant_grid,... (1 : x_rawlength) * d_geant_grid, (1 : z_rawlength) * d_geant_grid); [y_kgrid, x_kgrid, z_kgrid] = meshgrid((1 : y_newlength) * kgrid.dy,... (1 : x_newlength) * kgrid.dx, (1 : z_newlength) * kgrid.dz); raw_source.p0 = interp3(y_geant, x_geant, z_geant, raw_source.p0,... y_kgrid, x_kgrid, z_kgrid); %% Put the source into the kgrid % This is similar to my makeSource_spatial3D function. It creates an empty % container of the size of the k-Wave grid, and puts the raw source into % this container. source.p0 = zeros(kgrid.Nx, kgrid.Ny, kgrid.Nz); y_shift = floor((kgrid.Ny - y_newlength) / 2); z_shift = floor((kgrid.Nz - z_newlength) / 2); source.p0(1 : x_newlength, y_shift + (1 : y_newlength),... z_shift + (1 : z_newlength)) = raw_source.p0; %% Creating sensor y_bragg = floor(kgrid.Ny/2)+1; % index z_bragg = floor(kgrid.Nz/2)+1; % index [~,x_bragg] = max(source.p0(:,y_bragg,z_bragg)); %index sensor.mask = zeros(Nx, Ny, Nz); sensor.mask(x_bragg + round(0.025/kgrid.dx), y_bragg, z_bragg) = 1; %% Additional input parameters input_fullpath = strcat('/project/med4/IONOAKUSTIK/SimulationAndreas/',... SaveToDisk,'.h5'); cvsmooth = [smooth_p0, smooth_soundspeed, smooth_density]; input_args={'PMLInside', PMLInside, 'Smooth', cvsmooth, 'DataCast',... DataCast, 'StreamToDisk', StreamToDisk, 'SaveToDisk', input_fullpath}; kspaceFirstOrder3D(kgrid,medium,source,sensor,input_args{:});</code></pre> Bradley Treeby on "Unexpected instability" http://www.k-wave.org/forum/topic/unexpected-instability#post-6517 Fri, 13 Jul 2018 21:54:49 +0000 Bradley Treeby 6517@http://www.k-wave.org/forum/ Hi Foufour, Can you create an example that doesn't depend on external functions or data that has the same problem? Then I can take a look. Brad. Foufour on "Unexpected instability" http://www.k-wave.org/forum/topic/unexpected-instability#post-6507 Thu, 28 Jun 2018 10:29:00 +0000 Foufour 6507@http://www.k-wave.org/forum/ Here comes the matlab code... <pre><code>%% Input parameters grid_length = 0.1; % length of kgrid in m grid_width = 0.1; % width of kgrid in m axial_spacing = 250e-6; % spacing between neighbouring points in axial direction (dx) in m lateral_spacing = 250e-6; % spacing between neighbouring points in lateral direction (dy,dz) in m time_step = 25e-9; % simulation time step in s total_time = 35e-6; % total simulation time in s T = 25; % temperature in °C source_filename = '/project/med4/IONOAKUSTIK/SimulationAndreas/Marco_SIRMIO_Geant4data/TestSignalForIonoacoustics_50MeV_SigmaEnergy_1MeV_voxelSize_0.1x0.1x0.1mm_FWHM_2mm.txt'; geant_grid_spacing = 100e-6; % grid spacing of the above source file detector_shape = 'point'; % String. Possible values: 'planar', 'spherical', 'point', 'rahmen', 'multiple_points', 'triple_planar' smooth_p0 = true; % Boolean controlling whether initial pressure distribution gets smoothed automatically smooth_density = true; % Boolean controlling whether density distribution gets smoothed automatically smooth_soundspeed = true; % Boolean controlling whether sound speed distribution gets smoothed automatically StreamToDisk = true; % Boolean controlling whether simulation data are stored on disk during simulation DataCast = 'gpuArray-double'; % Variable precision & simulation mode. Possible: 'single', 'double', 'gpuArray-double', 'gpuArray-single' SaveToDisk = 'sirmiobeam-double-250-in'; % File name of the HDF5 input file PMLInside = false; %% Creating kWaveGrid % In all dimensions, the grid should have as small prime factors as % possible. That's why I chose to increase the grid to the next power of 2. % The 20 points subtracted will get re-added in the end as PML. Nx = 2^nextpow2(grid_length/axial_spacing)-20; Ny = 2^nextpow2(grid_width/lateral_spacing)-20; Nz = Ny; kgrid = kWaveGrid(Nx,axial_spacing,Ny,lateral_spacing,Nz,lateral_spacing); kgrid.dt = time_step; kgrid.Nt = round(total_time/time_step); %% Creating medium sound_speedwater=1.40238677e3 + 5.03798765*T -5.80980033e-2*T.^2 + 3.34296650e-4*T.^3 -1.47936902e-6*T.^4 + 3.14893508e-9*T.^5; %[m/s] densitywater = 999.9+T.*(0.06079+T.*(-0.008313+T.*(6.595e-5+T.*(-4.437e-7+1.369e-9.*T))));%[kg/m^3] medium = struct('sound_speed', sound_speedwater, 'density', densitywater, 'sound_speed_ref', sound_speedwater); %% Crating source source = new_Read_Geant(source_filename,T); % Get the data from *.dat format into matlab [source, ~] = adjust_grids(kgrid,source,geant_grid_spacing,0); % I added this because I wanted to become independent of the grid spacing % of the Geant4 files. So I decided to transform the source file from its % original grid spacing into that of my kgrid source = makeSource_spatial3D(source,kgrid.Nx,kgrid.Ny,kgrid.Nz,kgrid.dx,0,0,0,0,0); % Put the source into the kgrid. %% Creating sensor y_bragg = floor(kgrid.Ny/2)+1; % index z_bragg = floor(kgrid.Nz/2)+1; % index [~,x_bragg] = max(source.p0(:,y_bragg,z_bragg)); %index sensor.mask = zeros(Nx, Ny, Nz); sensor.mask(x_bragg + round(0.025/kgrid.dx), y_bragg, z_bragg) = 1; %% Additional input parameters input_fullpath = strcat('/project/med4/IONOAKUSTIK/SimulationAndreas/',SaveToDisk,'.h5'); cvsmooth = [smooth_p0, smooth_soundspeed, smooth_density]; input_args={'PMLInside', PMLInside, 'Smooth', cvsmooth, 'DataCast', DataCast, 'StreamToDisk', StreamToDisk, 'SaveToDisk', input_fullpath}; kspaceFirstOrder3D(kgrid,medium,source,sensor,input_args{:});</code></pre> Foufour on "Unexpected instability" http://www.k-wave.org/forum/topic/unexpected-instability#post-6506 Thu, 28 Jun 2018 09:45:10 +0000 Foufour 6506@http://www.k-wave.org/forum/ Hi, I'm facing instability issues I don't understand. I'm simulating in 3D in homogeneous water. On coarse grids (1 mm and 500 µm spacings, with time steps of 100 and 50 ns, respectively), everything behaves exactly as I expect. However, when I set the spacing to 250 µm, the simulation is no longer stable, but blows up immediately, no matter how small I choose the time step (even tried 1 ns). The medium doesn't change (since it's homogeneous), so I checked the source. I'm using the same source file for all grids and match it onto the kgrid. This happens automatically, the only code lines I change are really the grid spacing and the time step. I re-imported the source from the HDF5 input file to investigate it. But it looks fine. No steeper gradients than on the coarse grid (max. 25 Pa/m, corresponding to 6 mPa between two grid points), and the shape looks fine as well. I'm running out of explanations... I'll extract the matlab code lines for you. Could you please give me an idea of what went wrong? Kind regards, Foufour