k-Wave User Forum » User Favorites: Elseliv

k-Wave User Forum » User Favorites: Elseliv elseliv Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 03:48:29 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php HD Y on "Grid size and magnitude of pressure" http://www.k-wave.org/forum/topic/grid-size-and-magnitude-of-pressure#post-9187 Wed, 19 Feb 2025 23:18:09 +0000 HD Y 9187@http://www.k-wave.org/forum/ Bradley Treeby on "Validation of the tissue velocity recorded by sensor" http://www.k-wave.org/forum/topic/validation-of-the-tissue-velocity-recorded-by-sensor#post-8125 Mon, 19 Apr 2021 13:46:57 +0000 Bradley Treeby 8125@http://www.k-wave.org/forum/ Hi Ekaterina, I couldn't view the image, however, if you experience an instability you could try reducing the time step size. On the topic of reflections, I'm not sure these will match any physically realistic scenario, as you are imposing a time-varying boundary condition, thus the reflections will depend on the values for the velocity that you are imposing at each time step in the simulation. Brad. Elseliv on "Validation of the tissue velocity recorded by sensor" http://www.k-wave.org/forum/topic/validation-of-the-tissue-velocity-recorded-by-sensor#post-8108 Thu, 25 Mar 2021 19:13:33 +0000 Elseliv 8108@http://www.k-wave.org/forum/ I actually tried using the Dirichlet condition but this brought instability and the simulation got totally crashed. Here is a link with 2 plots for Dirichlet and NoDirichlet cases <a href="https://we.tl/t-Dft5lKj5vd" rel="nofollow">https://we.tl/t-Dft5lKj5vd</a> and it depicts in the top row - source velocities as I imposed them and on the 2-nd and 3-d rows - velocities measure by the sensor exactly next to the source. Could you then assume what might be a reason for such a crash? Concerning reflections, the source I use is located at the boundary between tissue and water so reflections are expected to happen anyway. Bradley Treeby on "Validation of the tissue velocity recorded by sensor" http://www.k-wave.org/forum/topic/validation-of-the-tissue-velocity-recorded-by-sensor#post-8106 Thu, 25 Mar 2021 17:52:22 +0000 Bradley Treeby 8106@http://www.k-wave.org/forum/ Hi Ekaterina, By default, a velocity source is a force source. It will inject or "add" to the field that is already there, and thus the actual velocity will be higher than the velocity you specify. If you want the source to act like a boundary condition (exactly impose the values you specify), then you can set <code>source.u_mode = 'dirichlet'</code>. Keep in mind that in this case, the imposed boundary condition will also reflect any ways coming back to the source. Brad. Elseliv on "Validation of the tissue velocity recorded by sensor" http://www.k-wave.org/forum/topic/validation-of-the-tissue-velocity-recorded-by-sensor#post-8091 Wed, 24 Mar 2021 01:32:51 +0000 Elseliv 8091@http://www.k-wave.org/forum/ Hello Dr.Treeby, Thank you for the great toolbox! I`m doing shear wave simulation in a thin tissue layer using pstdElastic2D. As a source.ux, I use a line where each element imposes one cycle of a sine wave with a frequency of 50 Hz and max amplitude of 0.03 m/s. Then I record u_split_field with a line of sensors perpendicular to the source line and get values significantly higher (sensor_data.ux_split_s ~ 0.1 m/s sensor_data.uy_split_s ~ 0.5 m/s), than those imposed by source. Even for the sensors placed tight to the source. What could become a reason for differences in the velocities imposed and observed? Is there a way to validate the magnitude of output data? May there be a problem with output units or something else? Best regards, Ekaterina Bradley Treeby on "Want to run a program on a GPU" http://www.k-wave.org/forum/topic/want-to-run-a-program-on-a-gpu#post-7449 Fri, 01 May 2020 08:14:35 +0000 Bradley Treeby 7449@http://www.k-wave.org/forum/ Hi Jack, Open the MATLAB help browser and have a look at the examples under the "Using the C++ code" heading. Brad. jackYANG on "Want to run a program on a GPU" http://www.k-wave.org/forum/topic/want-to-run-a-program-on-a-gpu#post-7446 Thu, 30 Apr 2020 04:04:26 +0000 jackYANG 7446@http://www.k-wave.org/forum/ Hello everyone, since my simulation model is relatively large, it takes a long time to run on the CPU, now I want to switch to running on the GPU.Can someone tell me in detail how my kspaceFirstOrder3D program works on a GPU. rainbowwhale on "Grid size and magnitude of pressure" http://www.k-wave.org/forum/topic/grid-size-and-magnitude-of-pressure#post-7389 Mon, 13 Apr 2020 10:25:44 +0000 rainbowwhale 7389@http://www.k-wave.org/forum/ Thank you for kind answer. I may confused between Pressure and power. I'll tried to integrate the power and make comparisons. bencox on "Grid size and magnitude of pressure" http://www.k-wave.org/forum/topic/grid-size-and-magnitude-of-pressure#post-7381 Fri, 10 Apr 2020 10:52:42 +0000 bencox 7381@http://www.k-wave.org/forum/ Hi rainbowwhale, Yes, you will see a change in the amplitude for a couple of reasons. The most obvious is that you are inputting the source on more grid points per unit length so the amplitude will be correspondingly higher. (Think of the sources as point sources on a collection of grid points, not sources filling voxels ie. not spread continuously over an area.) The second is that the grid can support higher frequencies as you reduce the grid spacing, and as your <code>source.p</code> signal contains a broad range of frequencies, more of the higher frequencies will be propagated in the grids with smaller grid spacing<code>dG</code>. Final point. You're currently changing the size of the sensor and the source simultaneously. It might be worth investigating the size of the source with the sensor as a single grid point, and then subsequently to investigate the effect of having a larger sensor area. Hope that helps, Ben rainbowwhale on "Grid size and magnitude of pressure" http://www.k-wave.org/forum/topic/grid-size-and-magnitude-of-pressure#post-7340 Wed, 01 Apr 2020 07:30:03 +0000 rainbowwhale 7340@http://www.k-wave.org/forum/ Hi I'm new to k-Wave and have a question about grid size I made a source shaped in square with same size for different grid size. I thought all cases have same area of source and result in same pressure magnitude but small grid case which have more grid cells in source showed higher pressure. Magnitude of pressure is proportional to the number of grid per length (not area!) and thus increased linearly with respect to sqrt(number of grid in source area). I have tried to use 'dirichlet' type source but showed same results. --> There was an typo and with a correct command, 'dirichlet' type worked. I'll appreciate your kind advises. Thank you %% clearvars clc PML_Size = 20; PML_Alpha = 2; input_args = {'PMLSize', PML_Size, 'PMLInside', true, 'PMLAlpha', PML_Alpha}; %% CFL = 0.15; dG = 0.01*1e-3; % reference grid size dt = CFL*dG/1500; Nt = 5000; Nx = 512; Nz = 512; pulse_Rect = conv(ones(1, ceil(40e-9/dt)), ones(1,7)/7); Tx_Rect = [pulse_Rect, zeros(1, Nt - length(pulse_Rect))]; medium.sound_speed = 1500; medium.density = 1000; % number of grid per source (and sensor) NG = [2, 4, 8, 16, 32, 64]; results = zeros(Nt, length(NG)); for k=1:length(NG) dG = 0.16*1e-3/NG(k); kgrid = kWaveGrid(Nx, dG, Nz, dG); kgrid.dt = dt; kgrid.Nt = Nt; kgrid.t_array = (0:kgrid.Nt-1)*dt; %source.p_mode = 'dirichlet'; source.p_mask = zeros(Nx, Nz); source.p_mask((1:NG(k)) - 0.5*NG(k) + Nx/2, (1:NG(k)) - 0.5*NG(k) + Nz/2 - 3*NG(k)) = 1; source.p = Tx_Rect; sensor.mask = zeros(Nx, Nz); sensor.mask((1:NG(k)) - 0.5*NG(k) + Nx/2, (1:NG(k)) - 0.5*NG(k) + Nz/2 + 3*NG(k)) = 1; subplot(2, length(NG)+1, k) imagesc(kgrid.x_vec, kgrid.y_vec, (source.p_mask-sensor.mask)) axis image xlim([-1, 1]*0.9*1e-3) ylim([-1, 1]*0.9*1e-3) subplot(2, length(NG)+1, k + length(NG)+1) imagesc(kgrid.x_vec, kgrid.y_vec, (source.p_mask-sensor.mask)) axis image drawnow data = kspaceFirstOrder2DG(kgrid, medium, source, sensor, input_args{:}); results(:,k) = mean(data, 1)'; end subplot(2, length(NG)+1, length(NG) + 1) plot(results) Bradley Treeby on "Dirichlet condition doesnt work as expected" http://www.k-wave.org/forum/topic/dirichlet-condition-doesnt-work-as-expected#post-7212 Tue, 11 Feb 2020 12:00:44 +0000 Bradley Treeby 7212@http://www.k-wave.org/forum/ Hi grabka, I haven't quite managed to get to the bottom of your problem yet. First, I don't think you're plotting the same point. So I changed <code>sensor.mask = source.u_mask;</code> and then plotted <code>plot(kgrid.t_array, sensor_data.uz(1,:), kgrid.t_array, source.uz(1,:))</code> If you set the shear sound speed to zero, everything works as expected. However, when the shear speed is non-zero, the input and output don't agree. I'll keep digging and come back to you. Brad. grabka on "Dirichlet condition doesnt work as expected" http://www.k-wave.org/forum/topic/dirichlet-condition-doesnt-work-as-expected#post-7133 Thu, 21 Nov 2019 23:31:11 +0000 grabka 7133@http://www.k-wave.org/forum/ Hello. I have set up a simulation with homogeneous viscoelastic material and excite by a Velocity Gauss impulse in the middle. I have noticed that the velocities at the source positions are not equal to the defined ones, although I have defined dirichlet conditions. Why is that? In the manual it states: ‘dirichlet’: at each time step, the input pressure and velocity values are used to replace the existing values Here my code: <pre><code>%% variables f_max = 300; %[Hz] %% create the computational grid spacing = 5e-3; x_size = 0.15; %[m] total size x-direction y_size = 0.15; %[m] total size y-direction z_size = 0.014; %[m] total size z-direction dx = spacing; % grid point spacing in the x direction [m] dy = spacing; % grid point spacing in the y direction [m] dz = spacing; % grid point spacing in the z direction [m] Nx = ceil(x_size/dx)+2; % number of grid points in the x direction Ny = ceil(y_size/dy)+2; % number of grid points in the y direction Nz = ceil(z_size/dz)+2 % number of grid points in the z direction disp('Number of grid points: '); disp(Nx*Ny*Nz); kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); %% simulation time and step size dt = 0.1 * spacing / 1440.3; % CFL*dx/c_max CFL 0.3 default t_end = 0.5/300; % [s] kgrid.t_array = 0:dt:t_end; %% define the properties of the propagation medium medium.sound_speed_compression = ones(Nx,Ny,Nz) * 1440.3; % [m/s] medium.sound_speed_shear = ones(Nx,Ny,Nz) * 201.74; % [m/s] medium.alpha_coeff_compression = ones(Nx,Ny,Nz) * 0.5; % [10–3 m–1] medium.alpha_coeff_shear = ones(Nx,Ny,Nz) * 0.5; % [10–3 m–1] medium.density = ones(Nx,Ny,Nz) * 923.47; % [kg/m^3] %% define the source % define a square source element source_radius = 0.005/dx; % [grid points] source.u_mask = zeros(Nx,Ny,Nz); source.u_mask(Nx/2 - source_radius:Nx/2 + source_radius, Ny/2 - source_radius:Ny/2 + source_radius, Nz) = 1; nbr_of_sources = sum(source.u_mask,"all"); % define a gaussian impulse pos = gaussmf(kgrid.t_array, [t_end/20 t_end/6]); vel_gaussian = [0 diff(pos)/dt]; source.uz = repelem(vel_gaussian,nbr_of_sources, 1); source.u_mode= 'dirichlet'; %% define sensor sensor.mask = ones(Nx, Ny, Nz); %% run the simulation sensor.record = {'u'} % u particle velocity PML_size=5; % size of the PML in grid points input_args = {'PMLSize', PML_size, 'PlotPML', true, 'PMLInside', false}; sensor_data = pstdElastic3D(kgrid, medium, source, sensor, input_args{:});</code></pre> plotting the input vs ouput at sensor position gives different plots: plot(kgrid.t_array, sensor_data.uz(4659,:), kgrid.t_array, source.uz(1,:)) genny_p on "Instability in shear wave absorption simulation" http://www.k-wave.org/forum/topic/instability-in-shear-wave-absorption-simulation#post-6939 Fri, 28 Jun 2019 11:07:46 +0000 genny_p 6939@http://www.k-wave.org/forum/ Hi Brad, Thanks for the suggestion. The simulation seems to actually only be stable with longer time steps (sparser grid spacing and higher CFL number) and when I decrease the time step then the simulation becomes unstable (you can see this by changing the factors in the sample code I provided). I take it then this is not a result of a numerical instability. What other types of instabilities could be possible? Thanks, Genny xurui on "When will be released C++/CUDA GPU Simulation Codes?" http://www.k-wave.org/forum/topic/when-will-be-released-ccuda-gpu-simulation-codes#post-6934 Thu, 27 Jun 2019 14:59:49 +0000 xurui 6934@http://www.k-wave.org/forum/ Hi Brad, Thanks, I'm looking forward to it. I'll content myself with the Matlab version for now. Cheers, Rui. Bradley Treeby on "When will be released C++/CUDA GPU Simulation Codes?" http://www.k-wave.org/forum/topic/when-will-be-released-ccuda-gpu-simulation-codes#post-6900 Sat, 22 Jun 2019 20:56:27 +0000 Bradley Treeby 6900@http://www.k-wave.org/forum/ Hi Rui, This code is in a relatively stable beta, but won't be part of the upcoming V1.3 release. We hope to have it out later this year. Brad. Bradley Treeby on "Instability in shear wave absorption simulation" http://www.k-wave.org/forum/topic/instability-in-shear-wave-absorption-simulation#post-6894 Sat, 22 Jun 2019 20:18:49 +0000 Bradley Treeby 6894@http://www.k-wave.org/forum/ Hi Genny, For the shear code, we don't have an out of the box function to test for stability with an absorbing medium. For the fluid code, we have an approximate way of checking for the stable time step size (see the <code>checkStability</code> function). It might be possible to find a relationship for the elastic code, but I haven't looked at it. If it is a numerical instability, then keep making the CFL smaller, it should eventually work. Let me know if it doesn't. Brad. xurui on "When will be released C++/CUDA GPU Simulation Codes?" http://www.k-wave.org/forum/topic/when-will-be-released-ccuda-gpu-simulation-codes#post-6864 Wed, 29 May 2019 16:03:43 +0000 xurui 6864@http://www.k-wave.org/forum/ Hi, Will the GPU-accelerated 3D elastic wave propagation code (presented last year at the 2018 International Conference on High Performance Computing & Simulation) be released to the masses anytime soon? Cheers, Rui. genny_p on "Instability in shear wave absorption simulation" http://www.k-wave.org/forum/topic/instability-in-shear-wave-absorption-simulation#post-6851 Wed, 15 May 2019 12:05:54 +0000 genny_p 6851@http://www.k-wave.org/forum/ Hi Brad, Ben, and the k-Wave community, As I mentioned at Photonics West, I am applying k-Wave to model scattering of planar shear waves in heterogeneous tissue, and so far the toolbox has been very easy to use and is working quite well for elastic scattering problems. I recently started trying to model shear wave absorption using the medium.alpha_coeff_shear factor to account for viscoelasticity. The results look quite reasonable for certain cases, but when I increase the absorption factor enough, the result seems to become unstable. I am using a value on the order of 3e6 dB/(MHz^2 cm), since this is around what we measure for shear wave absorption in tissue-mimicking phantoms with wave frequency on the order of 600 Hz, i.e. 1.08 dB/cm. You can see the instability forming in the example code below, where the edges of the plane wave start to increase with an amplitude much larger than the initial wave amplitude. I saw in the Modelling Power Law Absorption Example that you discuss numerical errors (<a href="http://www.k-wave.org/documentation/example_na_modelling_absorption.php#heading3)" rel="nofollow">http://www.k-wave.org/documentation/example_na_modelling_absorption.php#heading3)</a>, but there it is mentioned that the errors can be minimized by reducing the time step. I found the instability I see seems to worsen with a decreased time step (which I implemented by reducing the CFL number). The instability is also worse when I reduce the grid spacing (which in turn also reduced the time step). See the comments in my code below for the factors that influence the instability that I investigated. All the parameters used in the example code are typical for our experiments, except maybe the inital velocity u_0, but there doesn't seem to be any dependence of the instability on this value. This leads me to the following questions: 1. Is this instability related to the numerical errors discussed in the Modelling Power Law Absorption Example? 2. Is there some relationship that tells me for which conditions (up to which absorption) the code works for before becoming unstable? 3. If I have a result that appears stable (i.e. with a low enough medium.alpha_coeff_shear value), how can I check whether the results are reliable, since when I decrease the grid size to check convergence the simulation then blows up? I would be grateful if someone had any insight to any of these questions. Thanks, Genny % ------------------------------------------------------------------------------------- % Modified from Shear Waves And Critical Angle Reflection Example % Simplified code that shows instability during absorption of planar shear wave % Heterogeneities removed to focus on instability % Grid size cropped for quicker simulation % The parameters coded below result in an amplitude that appears to blow up % May 2019 % % If one of the factors is changed to as decribed below (with all other parameters constant): % * decreasing alpha0_s1 to 2e6 the solution appears stable % * increasing alpha0_s1 to 4e6 the solution blows up more % * increasing t_end to 5000 the solution blows up more % * decreasing t_end to 2000 the solution appears stable % * decreasing cfl to 0.1 (decreases dt by a factor of 2) the solution blows up more % * increasing scale to 1.25 (decreases dx & dt by a factor of 0.8) the solution blows up much more % * decreasing scale to 0.75 (increases dx & dt by a factor of 1.3) the solution appears stable % * increasing cp1 to 3000 (decreases dt by a factor of 2) the solution blows up more % * decreasing cp1 to 750 (increases dt by a factor of 2) the solution appears stable % * increasing cs1 to 4.5 the solution blows up more % * decreasing cs1 to 3.5 the solution appears stable % * instability appears insensitive to frequency (300-1800 Hz tested) % * instability appears insensitive to u_0 (0.5 to 1000 tested) close all clearvars % define the medium properties for the outise material cp1 = 1500; % compressional wave speed [m/s] cs1 = 4.35; % shear wave speed [m/s] rho1 = 1000; % density [kg/m^3] alpha0_p1 = 0; % compressional absorption [dB/(MHz^2 cm)] alpha0_s1 = 3e6; % shear absorption [dB/(MHz^2 cm)] (measured approx 1.08 dB/cm at 600 Hz) % create the computational grid scale = 1; % grid spacing factor PML_size = 10; % [grid points] Nx = 40*scale ; % [grid points] Ny = 20*scale ; % [grid points] dx = 0.4e-3/scale; % [m] dy = 0.4e-3/scale; % [m] kgrid = kWaveGrid(Nx, dx, Ny, dy); % create the time array to run until wave passes through entire grid cfl = 0.2; t_end = 4500*1e-6; kgrid.makeTime(max([cp1 cs1]), cfl, t_end); % define the source source_freq = 600; % [Hz] u_0 = 500; % initial velocity source_mask = zeros(Nx, Ny); source_mask(:,1) = 1; source.u_mask = source_mask; source.ux = u_0 * sin(2 * pi * source_freq * kgrid.t_array); source.uy = 0 * sin(2 * pi * source_freq * kgrid.t_array); % define sensor sensor.record = {'u_max_all','u_min_all','u_final'}; % define the medium properties clear medium medium.sound_speed_compression = cp1*ones(Nx, Ny); medium.sound_speed_shear = cs1*ones(Nx, Ny); medium.density = rho1*ones(Nx, Ny); medium.alpha_coeff_compression = alpha0_p1*ones(Nx, Ny); medium.alpha_coeff_shear = alpha0_s1*ones(Nx, Ny); % set the input arguments input_args = {'PMLSize', PML_size, 'PMLAlpha', 2, 'PlotPML', false, ... 'PMLInside', false, 'PlotScale', 'auto', ... 'DataCast', 'single'}; % run the elastic simulation sensor_data_elastic = pstdElastic2D(kgrid, medium, source, sensor, input_args{:}); % visualization ui = u_0*cp1/cs1; % initial velocity figure(1) subplot(2,2,1); imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.ux_max_all/ui); caxis([0 1.5]) xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_x');title('ux max all') subplot(2,2,2); imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.ux_final/ui); caxis([-1.5 1.5]) xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_x');title('ux final') subplot(2,2,3); imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.uy_max_all/ui); caxis([0 1.5]) xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_y');title('uy max all') subplot(2,2,4); imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.uy_final/ui); caxis([-1.5 1.5]) xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_y');title('uy final') Alexll7 on "Multi-GPU fluid version and Elastic GPU version" http://www.k-wave.org/forum/topic/multi-gpu-fluid-version-and-elastic-gpu-version#post-6591 Thu, 04 Oct 2018 14:55:23 +0000 Alexll7 6591@http://www.k-wave.org/forum/ Hi Bradley, Thank you for your answer and the excellent work that you're doing for K_Wave. I hope these codes will come out soon. Best regards, Alex Bradley Treeby on "Multi-GPU fluid version and Elastic GPU version" http://www.k-wave.org/forum/topic/multi-gpu-fluid-version-and-elastic-gpu-version#post-6580 Thu, 13 Sep 2018 16:00:35 +0000 Bradley Treeby 6580@http://www.k-wave.org/forum/ Hi Alex, The multi-GPU fluid code is finished in regards to capability, but not finished in regards to usability. The current beta requires a reasonable number of command line variables, for example, to bind GPUs and threads to cores, to pad the subdomains to give small prime factors, and so on. Filip (the main programmer for this code) has started working on automating all of this, but it still might be some months before a release. The GPU elastic code is also working, however, Kristian (the main programmer for this code) is currently implementing higher-order time-stepping schemes to improve performance. Again, probably some months away before a release. Brad.