Commit 8e5963d0 authored by tilman.metz's avatar tilman.metz

added feature: rheinstetten model

parent 573b32d7
%3D plot of model
clear all;
clc
close all
fname='rheinstetten.';
parameter='vs';
source_file='../../par/sources/sources_rheinstetten.dat';
rec_file='../../par/receiver/receiver_rheinstetten.dat';
iteration=60;
nx=220; ny=100; nz=220; %ny:vertical
outx=1; outy=1; outz=1;
dh=0.15;
FW=10;
nx=nx/outx;ny=ny/outy;nz=nz/outz;
X=dh:dh*outx:nx*outx*dh;
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;
if strcmp(parameter,'vs')
caxis_value_1=0;%vs
caxis_value_2=300;%vs
end
if strcmp(parameter,'vp')
caxis_value_1=0;%vp
caxis_value_2=1700;%vp
end
% No changes to be done after this line---------
%------------------------------------------------
%Input of receiver and source positions--
fp = fopen(source_file);
SRC=textscan(fp, '%f %f %f %f %f %f', 'delimiter','\t', 'commentstyle','%','MultipleDelimsAsOne',1);
fclose(fp);
%
% fp = fopen(rec_file);
% REC=textscan(fp, '%f %f %f ', 'delimiter','\t', 'commentstyle','#','MultipleDelimsAsOne',1);
% fclose(fp);
% Lines to mark PML position
fontsize=14;
x_line=[FW FW nx-FW nx-FW FW]*dh;
y_line=[FW ny-FW ny-FW FW FW]*dh;
z_line=[FW nz-FW nz-FW FW FW]*dh;
%m --> model indices
slizex=int64(slx/dh/outx);
slizey=int64(sly/dh/outy);
slizez=int64(slz/dh/outz);
file_inp1=['../../par/model/' fname parameter '_it0'];
%file_inp2=['../../par/model/' fname parameter '_it' num2str(iteration)];
%Read 3D Models------------------------
fid=fopen(file_inp1,'r','ieee-le');
modelvec=zeros(ny/outy,nx/outx,nz/outz);
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);
%Creation of 2D Slices--------------
model_truex=zeros(ny,nz);
model_truey=zeros(nz,nx);
model_truez=zeros(ny,nx);
for y=1:1:ny
for z=1:1:nz
model_truex(y,z)=model_true(y,slizex,z);
end
end
for x=1:1:nx
for z=1:1:nz
model_truey(z,x)=model_true(slizey,x,z);
end
end
for x=1:1:nx
for y=1:1:ny
model_truez(y,x)=model_true(y,x,slizez);
end
end
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
%Plotting-------------------------------------------
figure(1)
imagesc(Z,Y,model_truex);
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('z in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['true ' parameter '-model. 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{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{1}(ii)==slx;
% plot(REC{3}(ii),REC{2}(ii),'.blue','MarkerSize',15,'LineWidth',2.0);
% end
%
% end
hold off
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');
%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{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{1}(ii)==slx;
% plot(REC{3}(ii),REC{2}(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')
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(['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{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{2}(ii)==sly;
% plot(REC{1}(ii),REC{3}(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')
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)
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);
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);
% end
%
% end
hold off
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('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{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{3}(ii)==sly;
% plot(REC{1}(ii),REC{2}(ii),'.blue','MarkerSize',15,'LineWidth',2.0);
% end
%
% end
hold off
figure(6)
imagesc(X,Y,modelz);
line(x_line,y_line,'LineStyle','--','Color','k')
colb=colorbar;
coll=get(colb,'xlabel');
colormap('jet');
set(coll,'String',[parameter ' in m/s'],'FontSize',fontsize);
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([parameter '-model. Iteration ' num2str(iteration) ' 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{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{3}(ii)==sly;
% plot(REC{1}(ii),REC{2}(ii),'.blue','MarkerSize',15,'LineWidth',2.0);
% end
%
% end
hold off
%exportfig(2, [fname parameter '.eps'],'bounds','tight', 'color','rgb', ...
% 'preview','none', 'resolution',200, 'lockaxes',1);
#FLNODE
#converted by MOCOX V1.2 from final.mod.gremlin
# 0.0000 1 0 6371.0000
# depth rho vp vs
# in m in g/cm^3 in km/s in km/s
#7
0.0000 1.7030 0.3569 0.0917
0.9299 1.7030 0.3569 0.1661
0.9299 1.6960 0.3569 0.1661
6.2636 1.6960 0.3569 0.2479
6.2636 2.0010 1.7450 0.2735
6.3800 2.0010 1.7450 0.2751
6.3800 2.0010 1.7450 0.2751
{
"MODELING PARAMETERS" : "comment",
"Note that y denotes the vertical direction !" : "comment",
"Domain Decomposition" : "comment",
"NPROCX" : "2",
"NPROCY" : "2",
"NPROCZ" : "2",
"3-D Grid" : "comment",
"NX" : "220",
"NY" : "100",
"NZ" : "220",
"DX" : "0.15",
"DY" : "0.15",
"DZ" : "0.15",
"FD order" : "comment",
"FDORDER" : "4",
"FDCOEFF" : "2",
"Time Stepping" : "comment",
"TIME" : "0.6",
"DT" : "4.0e-05",
"Source" : "comment",
"SOURCE_SHAPE" : "4",
"SOURCE_TYPE" : "4",
"SRCREC" : "1",
"SOURCE_FILE" : "./sources/sources_rheinstetten.dat",
"RUN_MULTIPLE_SHOTS" : "1",
"Model" : "comment",
"READMOD" : "0",
"MFILE" : "model/rheinstetten",
"Q-approximation" : "comment",
"L" : "0",
"Free Surface" : "comment",
"FREE_SURF" : "1",
"Absorbing Boundary" : "comment",
"ABS_TYPE" : "1",
"FW" : "10",
"VPPML" : "1700.0",
"FPML" : "31.25",
"Snapshots" : "comment",
"SNAP" : "0",
"Receiver" : "comment",
"SEISMO" : "1",
"READREC" : "2",
"Receiver array" : "comment",
"REC_ARRAY" : "1",
"REC_ARRAY_DEPTH" : "0.15",
"REC_ARRAY_DIST" : "30.0",
"DRX" : "10",
"DRZ" : "10",
"Seismograms" : "comment",
"NDT" : "1",
"SEIS_FORMAT" : "1",
"SEIS_FILE" : "./su_obs/obs_rheinstetten",
"Method" : "comment",
"METHOD" : "0",
"MOD_OUT_FILE" : "./model/rheinstetten"
}
# run forward simulation to obtain observed data
#make clean
make install MODEL=../genmod/rheinstetten_elastic.c
mpirun -np 8 nice -19 ../bin/ifos3d ./in_and_out/rheinstetten/rheinstetten_FW.json | tee ./in_and_out/rheinstetten/ifos3D.out
cp model/rheinstetten.vs_it0 model/rheinstetten.vs.true
cp model/rheinstetten.vp_it0 model/rheinstetten.vp.true
cp model/rheinstetten.rho_it0 model/rheinstetten.rho.true
# invert observed data using homogeneous starting model
#make clean
#make install MODEL=../genmod/rheinstetten_start.c
# mpirun -np 8 nice -19 ../bin/ifos3d ./in_and_out/rheinstetten/ifos3d_rheinstetten.json | tee ./in_and_out/rheinstetten/ifos3D_INV.out
3.75 0.15 3.75 0.000 31.5 1.0
3.75 0.15 8.25 0.000 31.5 1.0
3.75 0.15 12.75 0.000 31.5 1.0
3.75 0.15 17.25 0.000 31.5 1.0
3.75 0.15 21.75 0.000 31.5 1.0
3.75 0.15 26.25 0.000 31.5 1.0
8.25 0.15 3.75 0.000 31.5 1.0
8.25 0.15 8.25 0.000 31.5 1.0
8.25 0.15 12.75 0.000 31.5 1.0
8.25 0.15 17.25 0.000 31.5 1.0
8.25 0.15 21.75 0.000 31.5 1.0
8.25 0.15 26.25 0.000 31.5 1.0
3.75 0.15 3.75 0.000 31.5 1.0
12.75 0.15 8.25 0.000 31.5 1.0
12.75 0.15 12.75 0.000 31.5 1.0
12.75 0.15 17.25 0.000 31.5 1.0
12.75 0.15 21.75 0.000 31.5 1.0
12.75 0.15 26.25 0.000 31.5 1.0
17.25 0.15 3.75 0.000 31.5 1.0
17.25 0.15 8.25 0.000 31.5 1.0
17.25 0.15 12.75 0.000 31.5 1.0
17.25 0.15 17.25 0.000 31.5 1.0
17.25 0.15 21.75 0.000 31.5 1.0
17.25 0.15 26.25 0.000 31.5 1.0
21.75 0.15 3.75 0.000 31.5 1.0
21.75 0.15 8.25 0.000 31.5 1.0
21.75 0.15 12.75 0.000 31.5 1.0
21.75 0.15 17.25 0.000 31.5 1.0
21.75 0.15 21.75 0.000 31.5 1.0
21.75 0.15 26.25 0.000 31.5 1.0
26.25 0.15 3.75 0.000 31.5 1.0
26.25 0.15 8.25 0.000 31.5 1.0
26.25 0.15 12.75 0.000 31.5 1.0
26.25 0.15 17.25 0.000 31.5 1.0
26.25 0.15 21.75 0.000 31.5 1.0
26.25 0.15 26.25 0.000 31.5 1.0
%
% Definition of (distributed) source positions:
% position, time-delay, centre frequency, and amplitude
%
% (comment line is indicated by # as first character)
%
% Parameters for each source node (one source node per line):
%
% NSRC (one separateline)
% XSRC YSRC ZSRC TD FC AMP
%
% Symbols:
% XSRC= x-coordinate of source point [meter]
% YSRC= y-coordinate of source point [meter] (vertical)
% ZSRC= z-coordinate of source point [meter]
% TD= excitation time (time-delay) for source node [s]
% FC= centre frequency of source signal [Hz]
% AMP= maximum amplitude of source signal
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