Commit 4afe7b37 authored by tilman.metz's avatar tilman.metz

implemented trench in rheinstetten genmod file

parent fe5c2afc
......@@ -10,7 +10,7 @@ void model(st_model *mod){
/* extern variables */
extern int NX, NY,NZ, NXG, NYG, NZG, POS[4], MYID;
extern float DZ;
extern float DX,DY,DZ;
extern FILE *FP;
/* local variables */
......@@ -23,6 +23,39 @@ void model(st_model *mod){
float *fldepth, *flrho, *flvp, *flvs;
/* rench=1 - V shaped trench */
int trench=1;
float width=10;
float depth=3;
float pos=NXG/2*DX;
// Parameters of trench "Ettlinger Linie" velocities in km/s and density in kg/m^3
float evs=0.0917;
float evp=0.2569;
float erho=1.5;
/* 2 planes as border for the v-shape */
/* hess-form: n+(x-a)=0 (n=normal vektor a=point of plane)*/
float n1[4], a1[4],n2[4], a2[4], dist,rest1,rest2;
dist=sqrt(0.5*width*width);
a1[1]=dist;
a1[2]=0;
a1[3]=0;
a2[1]=0;
a2[2]=0;
a2[3]=dist;
n1[1]=-depth;
n1[2]=-dist;
n1[3]=depth;
n2[1]=-depth;
n2[2]=dist;
n2[3]=depth;
nodes=7;
fldepth=vector(1,nodes);
......@@ -73,15 +106,15 @@ void model(st_model *mod){
if ((i==1) && (k==1) && (MYID==0)){
fprintf(FP,"depth: %f m: double node\n",fldepth[l]);}}
else{
for (j=(int)(fldepth[l]/DZ)+1;j<=(int)(fldepth[l+1]/DZ);j++){
for (j=(int)(fldepth[l]/DY)+1;j<=(int)(fldepth[l+1]/DY);j++){
vp=0.0; vs=0.0; rhov=0.0;
vp=(DZ*(j-1)-fldepth[l])*(flvp[l+1]-flvp[l])/(fldepth[l+1]-fldepth[l])+flvp[l];
vp=(DY*(j-1)-fldepth[l])*(flvp[l+1]-flvp[l])/(fldepth[l+1]-fldepth[l])+flvp[l];
vp=vp*1000.0;
vs=(DZ*(j-1)-fldepth[l])*(flvs[l+1]-flvs[l])/(fldepth[l+1]-fldepth[l])+flvs[l];
vs=(DY*(j-1)-fldepth[l])*(flvs[l+1]-flvs[l])/(fldepth[l+1]-fldepth[l])+flvs[l];
vs=vs*1000.0;
rhov=(DZ*(j-1)-fldepth[l])*(flrho[l+1]-flrho[l])/(fldepth[l+1]-fldepth[l])+flrho[l];
rhov=(DY*(j-1)-fldepth[l])*(flrho[l+1]-flrho[l])/(fldepth[l+1]-fldepth[l])+flrho[l];
rhov=rhov*1000.0;
......@@ -106,7 +139,7 @@ void model(st_model *mod){
}
}
for (j=(int)(fldepth[nodes]/DZ)+1;j<=NYG;j++){
for (j=(int)(fldepth[nodes]/DY)+1;j<=NYG;j++){
vp=0.0; vs=0.0; rhov=0.0;
vp=flvp[nodes]*1000.0; vs=flvs[nodes]*1000.0; rhov=flrho[nodes]*1000.0;
......@@ -128,10 +161,45 @@ void model(st_model *mod){
mod->pi[jj][ii][kk]=piv;
}
}
/*Trench Ettlinger Linie*/
if (trench==1) {
for (j=1;j<=(int) depth/DY+1;j++){
rest1=n1[1]*(i*DX-a1[1])+n1[2]*(j*DY-a1[2])+n1[3]*(k*DZ-a1[3]);
rest2=n2[1]*(i*DX-a2[1])+n2[2]*(j*DY-a2[2])+n2[3]*(k*DZ-a2[3]);
if ((rest1>=0) && (rest2<=0)) {
vp=evp*1000.0;
vs=evs*1000.0;
rhov=erho*1000.0;
muv=vs*vs*rhov;
piv=vp*vp*rhov;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY)) &&
(POS[3]==((k-1)/NZ))) {
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=rhov;
mod->pi[jj][ii][kk]=piv;
}
}
}
}
}
}
free_vector(fldepth,1,nodes);
......
......@@ -9,6 +9,7 @@ parameter='vs';
source_file='../../par/sources/sources_rheinstetten.dat';
rec_file='../../par/receiver/receiver_rheinstetten.dat';
true_only=1;
iteration=60;
nx=220; ny=100; nz=220; %ny:vertical
......@@ -22,9 +23,9 @@ Y=dh:dh*outy:ny*outy*dh;
Z=dh:dh*outz:nz*outz*dh;
% Slizes of the 3D-Model in m
slx=10;
sly=10;
slz=10;
slx=5;
sly=2;
slz=5;
if strcmp(parameter,'vs')
caxis_value_1=0;%vs
......@@ -46,9 +47,11 @@ SRC=textscan(fp, '%f %f %f %f %f %f', 'delimiter','\t', 'commentstyle','%','Mult
fclose(fp);
%
% fp = fopen(rec_file);
% REC=textscan(fp, '%f %f %f ', 'delimiter','\t', 'commentstyle','#','MultipleDelimsAsOne',1);
% fclose(fp);
if (true_only==0)
fp = fopen(rec_file);
REC=textscan(fp, '%f %f %f ', 'delimiter','\t', 'commentstyle','#','MultipleDelimsAsOne',1);
fclose(fp);
end
% Lines to mark PML position
fontsize=14;
......@@ -64,8 +67,9 @@ slizez=int64(slz/dh/outz);
file_inp1=['../../par/model/' fname parameter '_it0'];
%file_inp2=['../../par/model/' fname parameter '_it' num2str(iteration)];
if (true_only==0)
file_inp2=['../../par/model/' fname parameter '_it' num2str(iteration)];
end
......@@ -77,11 +81,13 @@ modelvec=fread(fid,(nx*ny*nz),'float');
model_true=reshape(modelvec,ny,nx,nz);
% fid=fopen(file_inp2,'r','ieee-le');
% modelvec=zeros(ny/outy,nx/outx,nz/outz);
% modelvec=fread(fid,(nx*ny*nz),'float');
%
% model=reshape(modelvec,ny,nx,nz);
if (true_only==0)
fid=fopen(file_inp2,'r','ieee-le');
modelvec=zeros(ny/outy,nx/outx,nz/outz);
modelvec=fread(fid,(nx*ny*nz),'float');
end
model=reshape(modelvec,ny,nx,nz);
%Creation of 2D Slices--------------
......@@ -113,24 +119,25 @@ modelx=zeros(ny,nz);
modely=zeros(nz,nx);
modelz=zeros(ny,nx);
% for y=1:1:ny
% for z=1:1:nz
% modelx(y,z)=model(y,slizex,z);
% end
% end
%
% for x=1:1:nx
% for z=1:1:nz
% modely(z,x)=model(slizey,x,z);
% end
% end
%
% for x=1:1:nx
% for y=1:1:ny
% modelz(y,x)=model(y,x,slizez);
% end
% end
if (true_only==0)
for y=1:1:ny
for z=1:1:nz
modelx(y,z)=model(y,slizex,z);
end
end
for x=1:1:nx
for z=1:1:nz
modely(z,x)=model(slizey,x,z);
end
end
for x=1:1:nx
for y=1:1:ny
modelz(y,x)=model(y,x,slizez);
end
end
end
%Plotting-------------------------------------------
figure(1)
......@@ -164,124 +171,128 @@ end
% end
hold off
figure(2)
imagesc(Z,Y,modelx);
line(x_line,y_line,'LineStyle','--','Color','k')
figure(3)
imagesc(X,Z,model_truey);
line(x_line,z_line,'LineStyle','--','Color','k')
colb=colorbar;
coll=get(colb,'xlabel');
set(coll,'String',[parameter ' in m/s'],'FontSize',fontsize);
colormap('jet');
caxis([caxis_value_1 caxis_value_2])
set(gca,'ydir','reverse');
%set(gca,'xdir','reverse');
xlabel('z in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title([ parameter '-model. Iteration ' num2str(iteration) ' Slize at x=' num2str(slx) 'm'],'FontSize',fontsize)
set(gca,'ydir','normal');
xlabel('x in m','FontSize',fontsize)
ylabel('z in m','FontSize',fontsize)
title(['true ' parameter '-model. Slize at y=' num2str(sly) 'm'],'FontSize',fontsize)
axis tight
hold on
%Plot sources if in model slice
for ii=1:length(SRC{1})
if SRC{1}(ii)==slx;
plot(SRC{3}(ii),SRC{2}(ii),'*blue','MarkerSize',15,'LineWidth',2.0);
if SRC{2}(ii)==sly;
plot(SRC{1}(ii),SRC{3}(ii),'*blue','MarkerSize',15,'LineWidth',2.0);
end
end
%Plot receiver if in model slice
% for ii=1:length(REC{1})
% if REC{1}(ii)==slx;
% plot(REC{3}(ii),REC{2}(ii),'.blue','MarkerSize',15,'LineWidth',2.0);
% if REC{2}(ii)==sly;
% plot(REC{1}(ii),REC{3}(ii),'.blue','MarkerSize',15,'LineWidth',2.0);
% end
%
% end
hold off
figure(3)
imagesc(X,Z,model_truey);
line(x_line,z_line,'LineStyle','--','Color','k')
figure(5)
imagesc(X,Y,model_truez);
line(x_line,y_line,'LineStyle','--','Color','k')
colb=colorbar;
coll=get(colb,'xlabel');
set(coll,'String',[parameter ' in m/s'],'FontSize',fontsize);
colormap('jet');
caxis([caxis_value_1 caxis_value_2])
set(gca,'ydir','reverse');
%set(gca,'xdir','reverse');
xlabel('x in m','FontSize',fontsize)
ylabel('z in m','FontSize',fontsize)
title(['true ' parameter '-model. Slize at y=' num2str(sly) 'm'],'FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['true ' parameter '-model. Slize at z=' num2str(slz) 'm'],'FontSize',fontsize)
axis tight
hold on
%Plot sources if in model slice
for ii=1:length(SRC{1})
if SRC{2}(ii)==sly;
plot(SRC{1}(ii),SRC{3}(ii),'*blue','MarkerSize',15,'LineWidth',2.0);
if SRC{3}(ii)==slz;
plot(SRC{1}(ii),SRC{2}(ii),'*blue','MarkerSize',15,'LineWidth',2.0);
end
end
%Plot receiver if in model slice
% for ii=1:length(REC{1})
% if REC{2}(ii)==sly;
% plot(REC{1}(ii),REC{3}(ii),'.blue','MarkerSize',15,'LineWidth',2.0);
% if REC{3}(ii)==sly;
% plot(REC{1}(ii),REC{2}(ii),'.blue','MarkerSize',15,'LineWidth',2.0);
% end
%
% end
hold off
figure(4)
imagesc(X,Z,modely);
line(x_line,z_line,'LineStyle','--','Color','k')
if (true_only==0)
figure(2)
imagesc(Z,Y,modelx);
line(x_line,y_line,'LineStyle','--','Color','k')
colb=colorbar;
coll=get(colb,'xlabel');
set(coll,'String',[parameter ' in m/s'],'FontSize',fontsize);
colormap('jet');
caxis([caxis_value_1 caxis_value_2])
set(gca,'ydir','reverse');
xlabel('x in m','FontSize',fontsize)
ylabel('z in m','FontSize',fontsize)
title([parameter '-model. Iteration ' num2str(iteration) ' Slize at y=' num2str(sly) 'm'],'FontSize',fontsize)
%set(gca,'xdir','reverse');
xlabel('z in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title([ parameter '-model. Iteration ' num2str(iteration) ' Slize at x=' num2str(slx) 'm'],'FontSize',fontsize)
axis tight
hold on
%Plot sources if in model slice
for ii=1:length(SRC{1})
if SRC{2}(ii)==sly;
plot(SRC{1}(ii),SRC{3}(ii),'*blue','MarkerSize',15,'LineWidth',2.0);
if SRC{1}(ii)==slx;
plot(SRC{3}(ii),SRC{2}(ii),'*blue','MarkerSize',15,'LineWidth',2.0);
end
end
%Plot receiver if in model slice
% for ii=1:length(REC{1})
% if REC{2}(ii)==sly;
% plot(REC{1}(ii),REC{3}(ii),'.blue','MarkerSize',15,'LineWidth',2.0);
% if REC{1}(ii)==slx;
% plot(REC{3}(ii),REC{2}(ii),'.blue','MarkerSize',15,'LineWidth',2.0);
% end
%
% end
hold off
figure(5)
imagesc(X,Y,model_truez);
line(x_line,y_line,'LineStyle','--','Color','k')
figure(4)
imagesc(X,Z,modely);
line(x_line,z_line,'LineStyle','--','Color','k')
colb=colorbar;
coll=get(colb,'xlabel');
set(coll,'String',[parameter ' in m/s'],'FontSize',fontsize);
colormap('jet');
caxis([caxis_value_1 caxis_value_2])
set(gca,'ydir','reverse');
%set(gca,'xdir','reverse');
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['true ' parameter '-model. Slize at z=' num2str(slz) 'm'],'FontSize',fontsize)
ylabel('z in m','FontSize',fontsize)
title([parameter '-model. Iteration ' num2str(iteration) ' Slize at y=' num2str(sly) 'm'],'FontSize',fontsize)
axis tight
hold on
%Plot sources if in model slice
for ii=1:length(SRC{1})
if SRC{3}(ii)==slz;
plot(SRC{1}(ii),SRC{2}(ii),'*blue','MarkerSize',15,'LineWidth',2.0);
if SRC{2}(ii)==sly;
plot(SRC{1}(ii),SRC{3}(ii),'*blue','MarkerSize',15,'LineWidth',2.0);
end
end
%Plot receiver if in model slice
% for ii=1:length(REC{1})
% if REC{3}(ii)==sly;
% plot(REC{1}(ii),REC{2}(ii),'.blue','MarkerSize',15,'LineWidth',2.0);
% if REC{2}(ii)==sly;
% plot(REC{1}(ii),REC{3}(ii),'.blue','MarkerSize',15,'LineWidth',2.0);
% end
%
% end
......@@ -318,6 +329,11 @@ end
% end
hold off
end
%exportfig(2, [fname parameter '.eps'],'bounds','tight', 'color','rgb', ...
% 'preview','none', 'resolution',200, 'lockaxes',1);
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment