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);