http://www.k-wave.org/forum/topic/unstable-sim-when-reducing-grid-spacing-probably-my-fault Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 16:30:10 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/unstable-sim-when-reducing-grid-spacing-probably-my-fault#post-6136 Sun, 24 Sep 2017 11:04:01 +0000 Bradley Treeby 6136@http://www.k-wave.org/forum/ <p>Hi jwoodacr,</p> <p>When you change dx, you should also change dt. One way is to set this automatically using <code>makeTime</code>, e.g., </p> <p><code>kgrid.makeTime(medium.sound, CFL, t_end);</code></p> <p>Brad. </p> http://www.k-wave.org/forum/topic/unstable-sim-when-reducing-grid-spacing-probably-my-fault#post-6134 Fri, 22 Sep 2017 20:02:30 +0000 jwoodacr 6134@http://www.k-wave.org/forum/ <p>Hi,</p> <p>I've got a simple 1D simulation I'm running to examine ultrasound quarter-wavelength matching layers, and basically it works for some grid spacings (5um and above) and goes unstable for smaller ones (specifically 1um), which is a bit counter-intuitive to me. I am likely doing something wrong, but I can't seem to find how changing grid spacing would cause this.</p> <p>Note, the variable to change in the below code is dx. For a value of 5e-6 the sim works, and for 1e-6 it is unstable.</p> <pre><code>%% Create 1D Line Array dx = 5e-6; %grid point spacing freq = 5e6; %drive frequency of pressure pulse in Hz %Define Aluminum Properties alum.c = 6299; alum.d = 2698; alum.lambda = alum.c/freq; %Define Matching Layer Properties match_freq = 6.8e6; pary.c = 2135; pary.d = 1241.5; pary.lambda = pary.c/match_freq; %Define Water Properties water.c = 1540; water.d = 1000; water.lambda = pary.c/freq; %Define length of domains in wavelengths alum.num_wave = 10; water.num_wave = 10; %Calculate number of elements per domain, given dx and the above domain %lengths defined alum.n = floor(alum.num_wave*alum.lambda/dx); pary.n = floor(pary.lambda/(4*dx)); water.n = floor(water.num_wave*water.lambda/dx); Nx = alum.n+pary.n+water.n; %% Define Grid for 1D sim kgrid = makeGrid(alum.n+pary.n+water.n, dx); medium.density = [alum.d*ones(1,alum.n), pary.d*ones(1, pary.n), water.d*ones(1, water.n)]; medium.sound_speed = [alum.c*ones(1,alum.n), pary.c*ones(1, pary.n), water.c*ones(1, water.n)]; dt = 1e-9; %time step in seconds t_end = alum.num_wave+water.num_wave; %Time for wave to travel in both mediums t_end = t_end/freq; %Define time array t_array = 0:dt:t_end; kgrid.t_array = t_array; %Definine signal properties source_freq = freq; %Hz source_mag = 1; %Pa source_cycles = 4; %number of cycles %Define time array for source t_length = floor((source_cycles*(1/source_freq))/dt); t_length = min(t_length, length(kgrid.t_array)); %Define source values for pressure source source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array(1:t_length)); source.p = [source.p zeros(1, floor((2/source_freq)/dt))]; source.p = filterTimeSeries(kgrid, medium, source.p); %Define mask where pressure occurs source.p_mask = [zeros(1, floor(alum.n/2)),1]; source.p_mask = [source.p_mask, zeros(1, Nx-length(source.p_mask))]&#39;; sensor.mask = [floor(alum.n/2),alum.n+floor(water.n/2)].&#39;; sensor.record = {&#39;p_final&#39;, &#39;p_max&#39;, &#39;p_rms&#39;, &#39;p&#39;}; %% Run the Simultion display_mask = source.p_mask; Arraytype = &#39;single&#39;; %use gpuArray-single or gpuArray-double for gpu calcs %Arraytype = &#39;gpuArray-single&#39;; %use gpuArray-single or gpuArray-double for gpu calcs % assign the input options input_args = {&#39;DataCast&#39;, Arraytype,&#39;PlotSim&#39;, true, &#39;DisplayMask&#39;, display_mask, &#39;PMLInside&#39;, false, &#39;PlotPML&#39;, false}; %input_args = {&#39;DataCast&#39;, &#39;single&#39;,&#39;PlotSim&#39;, true, &#39;DisplayMask&#39;, display_mask, &#39;PMLInside&#39;, false, &#39;PlotPML&#39;, false, &#39;PlotScale&#39;, [-1 1]}; % run the simulation sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, input_args{:});</code></pre>