k-Wave User Forum » Topic: stress source lingers

k-Wave User Forum » Topic: stress source lingers http://www.k-wave.org/forum/topic/stress-source-lingers Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:29:16 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php ad000000 on "stress source lingers" http://www.k-wave.org/forum/topic/stress-source-lingers#post-6415 Mon, 16 Apr 2018 20:21:50 +0000 ad000000 6415@http://www.k-wave.org/forum/ Oops. Thanks. function myLine=myMakeLine(kgrid,lineStart,lineEnd); % function to make a line for K-wave, given cartesion input points. [lineStartGrid(1) lineStartGrid(2)] =cartPt2Grid(kgrid,lineStart); [lineEndGrid(1) lineEndGrid(2)]=cartPt2Grid(kgrid,lineEnd); myLine=makeLine(kgrid.Nx,kgrid.Ny,lineStartGrid,lineEndGrid); bencox on "stress source lingers" http://www.k-wave.org/forum/topic/stress-source-lingers#post-6411 Thu, 12 Apr 2018 20:38:43 +0000 bencox 6411@http://www.k-wave.org/forum/ Hi ad000000, I can't run your code as you have some homemade functions in there, but I recall seeing something similar in the past. Paste up your other functions and we'll take a look. Best wishes, Ben ad000000 on "stress source lingers" http://www.k-wave.org/forum/topic/stress-source-lingers#post-6387 Sat, 07 Apr 2018 02:48:49 +0000 ad000000 6387@http://www.k-wave.org/forum/ Hi. I made a linear stress source, to launch a plane wave at a hole to observe scattering. The source excitation is a wavelet. For some reason the simulation yields a finite stress (and velocity) along this line for the duration of the simulation. The BC's are normal (not Dirichlet). Any ideas? Thanks. CODE % HoleScatter.m % Launch a plane wave (shear of long) at a hole. clear source sensor medium kgrid %=== define variables===================== freq=5e6; epw=8; cp=5850; cs=3230; rho=8030; SDHdia=1.5e-3; type='long'; %shear or long dataInc=2; % Data is stored every dataInc grid points %==define geometry (remember, horizontal is Y direction)================== xlen=21e-3; ylen=16e-3; sdhXloc=0; sdhYloc=2e-3; srcYloc=-1e-3; %this shortens travel time, as opposed to putting it at left of grid. % ==================== setup grid ========================================= lambda=cp/freq; dx=lambda/epw; dy=lambda/epw; Nx= ceil(xlen/dx); Ny= ceil(ylen/dy); kgrid = makeGrid(Nx, dx, Ny, dy); %====set up material regions and kgrid ========================================= medium.sound_speed_compression =ones(kgrid.Nx, kgrid.Ny)*cp; medium.sound_speed_shear =ones(kgrid.Nx, kgrid.Ny)*cs; medium.alpha_coeff_compression = ones(kgrid.Nx, kgrid.Ny)*.3/(freq/1e6); %to give dB/(MHz^2 cm) medium.alpha_coeff_shear = ones(kgrid.Nx, kgrid.Ny)*.3/(freq/1e6)*4; medium.density=ones(kgrid.Nx, kgrid.Ny)*rho; SDHdiaGrid=SDHdia/dx; [SDHcenterGridx,SDHcenterGridy]=cartPt2Grid(kgrid,[sdhXloc sdhYloc]); disc = makeDisc(Nx, Ny, SDHcenterGridx, SDHcenterGridy, SDHdiaGrid/2); medium.sound_speed_compression(disc == 1) = 2000; medium.sound_speed_shear(disc == 1) = 1000; medium.density(disc==1)=rho;60; medium.alpha_coeff_compression(disc==1) = 1/(freq/1e6); %to give dB/(MHz^2 cm) medium.alpha_coeff_shear(disc==1) = 4/(freq/1e6)*4; if strcmp(type,'shear');vel=cs;else vel=cp;end t_end=(ylen/2+2*sdhYloc-srcYloc)/vel; %if ReverseFire==1; t_end=t_end/2;end dt=.3*dx/(cp); %set time steps (C should be Cmax) kgrid.t_array = 0:dt:t_end; Nt=length(kgrid.t_array); PML=20*dx; %======================= setup source =================================================== signal = wavelet(freq,4,0,2/freq,kgrid.t_array);%*1/(rho*cs); source.s_mask=myMakeLine(kgrid,[-.495*xlen, srcYloc],[.495*xlen, srcYloc ]); %source.s_mask=myMakeLine(kgrid,[-.495*xlen, srcYloc+3*SDHdia],[.495*xlen, srcYloc-3*SDHdia ]); srcinds=find(source.s_mask==1); if strcmp(type,'shear') source.sxy=repmat(signal,length(srcinds),1); else source.syy=repmat(signal,length(srcinds),1); end % The following is an attempt to use Dirichlet BC to simulate hole. Didn't work. %discinds=find(disc==1); %source.sxy=[source.sxy ;repmat(zeros(size(signal)),length(discinds),1)]; % source.u_mask=disc; % source.u_mask(:,1)=1; % srcinds=1:kgrid.Nx; % source.uy=repmat(signal,length(srcinds),1); % discinds=find(disc==1); % source.uy=[source.uy ;repmat(zeros(size(signal)),length(discinds),1)]; % source.u_mode = 'dirichlet'; %====setup sensor ===================================================== sensor.mask=zeros(kgrid.Nx, kgrid.Ny); sensor.mask(1:dataInc:kgrid.Nx,1:dataInc:kgrid.Ny)=1; sensor.record = {'u'}; %{'p','p_max'}; sensor.record_start_index=round((sdhYloc-srcYloc-SDHdia)/vel/dt); %===== display_mask ========================================= display_mask=disc; %==================== RUN ================================================= figpos=get(0,'DefaultFigurePosition'); set(0,'DefaultFigurePosition',[340,275,980,735]); input_args = {'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false,'PlotScale','auto','DataCast', 'single','RecordMovie',true,'MovieProfile','MPEG-4'}; sensor_data = pstdElastic2D(kgrid, medium, source, sensor, input_args{:}); %sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); set(0,'DefaultFigurePosition',figpos); kmodel=sensor_data; clear sensor_data medium source disc ; kmodel.freq=freq; kmodel.Nx=Nx;kmodel.Ny=Ny; kmodel.x_vec=kgrid.x_vec;kmodel.y_vec=kgrid.y_vec; kmodel.Nt=length(kgrid.t_array); kmodel.epw=epw; %uall=sqrt(kmodel.ux_max_all.^2+kmodel.uy_max_all.^2); u=sign(kmodel.uy).*sqrt(kmodel.ux.^2+kmodel.uy.^2); uh=abs(hilbert(u'))'; nox=ceil(kmodel.Nx/dataInc);noy=ceil(kmodel.Ny/dataInc); dat=zeros(nox,noy); for t=1:size(u,2); for j=1:nox for k=1:noy; dat(j,k,t)=uh((k-1)*nox+j,t); %dat(j,k,t)=u((k-1)*nox+j,t); %dat(j,k,t)=kmodel.ux((k-1)*nox+j,t); end end end % Record the movie %maxdat=max(max(max(dat))); %ax=[0 maxdat/500]; nframes=50; Nsamps=size(kmodel.ux,2); midframe=round(Nsamps/2); step=floor((Nsamps-1)/nframes); figure;imagesc(kmodel.y_vec,kmodel.x_vec,(dat(:,:,midframe)));axis equal;%caxis(ax) set(gca,'nextplot','replacechildren'); maxax=max(caxis); ax=[0 maxax/2]; for j= 1:nframes; t=1+(j-1)*step; imagesc(kmodel.y_vec,kmodel.x_vec,(dat(:,:,t)));caxis(ax) F(j) = getframe; end % Play the movie movie(F,2,4) %kmodel.uall=uall; kmodel.u=u; clear u save kmodel.mat kmodel -v7.3 %============================== function myLine=myMakeLine(kgrid,lineStart,lineEnd); % function to make a line for K-wave, given cartesion input points. [lineStartGrid(1) lineStartGrid(2)] =cartPt2Grid(kgrid,lineStart); [lineEndGrid(1) lineEndGrid(2)]=cartPt2Grid(kgrid,lineEnd); myLine=makeLine(kgrid.Nx,kgrid.Ny,lineStartGrid,lineEndGrid);