...
 
Commits (22)
......@@ -4,4 +4,5 @@
*.backup
*~*
*.orig
*.o
manual_IFOS3D.pdf
......@@ -24,8 +24,7 @@
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
void model(st_model *mod){
/*--------------------------------------------------------------------------*/
/* extern variables */
......@@ -55,7 +54,7 @@ float *** taus, float *** taup, float * eta){
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
mod->eta[l]=DT/pts[l];
}
ts=TAU;
tp=TAU;
......@@ -86,12 +85,12 @@ float *** taus, float *** taup, float * eta){
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
if(L){
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
}
......@@ -118,12 +117,12 @@ float *** taus, float *** taup, float * eta){
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
if(L){
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
......
......@@ -5,8 +5,7 @@
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float * eta){
void model(st_model *mod){
/*--------------------------------------------------------------------------*/
/* extern variables */
......@@ -36,7 +35,7 @@ float *** taus, float *** taup, float * eta){
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
mod->eta[l]=DT/pts[l];
}
ts=TAU;
tp=TAU;
......@@ -70,12 +69,13 @@ float *** taus, float *** taup, float * eta){
jj=j-POS[2]*NY;
kk=k-POS[3]*NZ;
if(L){
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
}
......@@ -92,7 +92,7 @@ float *** taus, float *** taup, float * eta){
sprintf(modfile,"%s.mod",MFILE);
writemod(modfile,rho,3);
writemod(modfile,mod,3);
MPI_Barrier(MPI_COMM_WORLD);
......
......@@ -5,8 +5,7 @@
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float *eta) {
void model(st_model *mod) {
/*-----------------------------------------------------------------------*/
/* extern variables */
......@@ -37,7 +36,7 @@ void model(float *** rho, float *** pi, float *** u,
for (l=1; l<=L; l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
mod->eta[l]=DT/pts[l];
}
ts=TAU;
......@@ -73,13 +72,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
}
......
......@@ -5,8 +5,7 @@
#include "fd.h"
void model(float *** rho, float *** pi, float *** u,
float *** taus, float *** taup, float *eta) {
void model(st_model *mod) {
/*-----------------------------------------------------------------------*/
/* extern variables */
......@@ -38,7 +37,7 @@ void model(float *** rho, float *** pi, float *** u,
for (l=1; l<=L; l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
mod->eta[l]=DT/pts[l];
}
ts=TAU;
......@@ -74,13 +73,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
}
......@@ -111,13 +110,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
......@@ -148,13 +147,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
......@@ -185,13 +184,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
......@@ -222,13 +221,13 @@ void model(float *** rho, float *** pi, float *** u,
kk=k-POS[3]*NZ;
if (L) {
taus[jj][ii][kk]=ts;
taup[jj][ii][kk]=tp;
mod->taus[jj][ii][kk]=ts;
mod->taup[jj][ii][kk]=tp;
}
u[jj][ii][kk]=muv;
rho[jj][ii][kk]=Rho;
pi[jj][ii][kk]=piv;
mod->u[jj][ii][kk]=muv;
mod->rho[jj][ii][kk]=Rho;
mod->pi[jj][ii][kk]=piv;
}
}
......
/*
* Model defined by flnode file.
*/
#include "fd.h"
void model(st_model *mod){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern int NX, NY,NZ, NXG, NYG, NZG, POS[4], MYID;
extern float DX,DY,DZ;
extern FILE *FP;
/* local variables */
float muv, piv, vp, vs, rhov;
int i, j,k, ii, jj,kk, l;
FILE *flfile;
int nodes;
char cline[256];
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.1000;
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);
flrho=vector(1,nodes);
flvp=vector(1,nodes);
flvs=vector(1,nodes);
/*read FL nodes from File*/
flfile=fopen("flnodes/flnodes.rheinstetten","r");
if (flfile==NULL) err(" FL-file could not be opened !");
for (l=1;l<=nodes;l++){
fgets(cline,255,flfile);
if (cline[0]!='#'){
sscanf(cline,"%f%f%f%f",&fldepth[l], &flrho[l], &flvp[l], &flvs[l]);
}
else l=l-1;
}
if(MYID==0){
fprintf(FP," ------------------------------------------------------------------ \n\n");
fprintf(FP," Information of FL nodes: \n\n");
fprintf(FP," \t depth \t rho \t vp \t vs \n\n");
for (l=1;l<=nodes;l++){
fprintf(FP," \t %f \t %f \t %f \t %f\n\n",fldepth[l],flrho[l],flvp[l],flvs[l]);
}
fprintf(FP," ------------------------------------------------------------------ \n\n");
}
/*-----------------------------------------------------------------------*/
/*-----------------------------------------------------------------------*/
/* loop over global grid */
for (i=1;i<=NXG;i++){
for (k=1; k<=NZG; k++) {
for (l=1;l<nodes;l++){
if (fldepth[l]==fldepth[l+1]){
if ((i==1) && (k==1) && (MYID==0)){
fprintf(FP,"depth: %f m: double node\n",fldepth[l]);}}
else{
for (j=(int)(fldepth[l]/DY)+1;j<=(int)(fldepth[l+1]/DY);j++){
vp=0.0; vs=0.0; rhov=0.0;
vp=(DY*(j-1)-fldepth[l])*(flvp[l+1]-flvp[l])/(fldepth[l+1]-fldepth[l])+flvp[l];
vp=vp*1000.0;
vs=(DY*(j-1)-fldepth[l])*(flvs[l+1]-flvs[l])/(fldepth[l+1]-fldepth[l])+flvs[l];
vs=vs*1000.0;
rhov=(DY*(j-1)-fldepth[l])*(flrho[l+1]-flrho[l])/(fldepth[l+1]-fldepth[l])+flrho[l];
rhov=rhov*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;
}
}
}
}
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;
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;
}
}
/*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);
free_vector(flrho,1,nodes);
free_vector(flvp,1,nodes);
free_vector(flvs,1,nodes);
}
/*
* Model defined by flnode file.
*/
#include "fd.h"
void model(st_model *mod){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern int NX, NY,NZ, NXG, NYG, NZG, POS[4], MYID;
extern float DZ;
extern FILE *FP;
/* local variables */
float muv, piv, vp, vs, rhov;
int i, j,k, ii, jj,kk, l;
FILE *flfile;
int nodes;
char cline[256];
float *fldepth, *flrho, *flvp, *flvs;
nodes=7;
fldepth=vector(1,nodes);
flrho=vector(1,nodes);
flvp=vector(1,nodes);
flvs=vector(1,nodes);
/*read FL nodes from File*/
flfile=fopen("flnodes/flnodes.rheinstetten_start","r");
if (flfile==NULL) err(" FL-file could not be opened !");
for (l=1;l<=nodes;l++){
fgets(cline,255,flfile);
if (cline[0]!='#'){
sscanf(cline,"%f%f%f%f",&fldepth[l], &flrho[l], &flvp[l], &flvs[l]);
}
else l=l-1;
}
if(MYID==0){
fprintf(FP," ------------------------------------------------------------------ \n\n");
fprintf(FP," Information of FL nodes: \n\n");
fprintf(FP," \t depth \t rho \t vp \t vs \n\n");
for (l=1;l<=nodes;l++){
fprintf(FP," \t %f \t %f \t %f \t %f\n\n",fldepth[l],flrho[l],flvp[l],flvs[l]);
}
fprintf(FP," ------------------------------------------------------------------ \n\n");
}
/*-----------------------------------------------------------------------*/
/*-----------------------------------------------------------------------*/
/* loop over global grid */
for (i=1;i<=NXG;i++){
for (k=1; k<=NZG; k++) {
for (l=1;l<nodes;l++){
if (fldepth[l]==fldepth[l+1]){
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++){
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=vp*1000.0;
vs=(DZ*(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=rhov*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;
}
}
}
}
for (j=(int)(fldepth[nodes]/DZ)+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;
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);
free_vector(flrho,1,nodes);
free_vector(flvp,1,nodes);
free_vector(flvs,1,nodes);
}
%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';
true_only=1;
iteration=60;
nx=220; ny=100; nz=220; %ny:vertical
outx=1; outy=1; outz=1;
dh=0.15;
FW=10;
Free_Surf=1;
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=2;
slz=20;
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);
%
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;
x_line=[FW FW nx-FW nx-FW FW]*dh;
z_line=[FW nz-FW nz-FW FW FW]*dh;
if(Free_Surf==1)
y_line=[0 ny-FW ny-FW 0 0]*dh;
else
y_line=[FW ny-FW ny-FW FW FW]*dh;
end
% Lines to mark slice position
x_line_slice=[slx slx];
y_line_slice=[sly sly];
z_line_slice=[slz slz];
%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'];
if (true_only==0)
file_inp2=['../../par/model/' fname parameter '_it' num2str(iteration)];
end
%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);
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--------------
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);
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)
imagesc(Z,Y,model_truex);
line(x_line,y_line,'LineStyle','--','Color','k')
line([0 nz*dh],y_line_slice,'LineStyle','-','Color','red')
%line(z_line_slice,[0 ny*dh],'LineStyle','-','Color','green')
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(3)
imagesc(X,Z,model_truey);
line(x_line,z_line,'LineStyle','--','Color','k')
line([0 nz*dh],z_line_slice,'LineStyle','-','Color','green')
%line(x_line_slice,[0 nx*dh],'LineStyle','-','Color','m')
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','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{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')
line([0 nz*dh],y_line_slice,'LineStyle','-','Color','red')
%line(x_line_slice,[0 ny*dh],'LineStyle','-','Color','magenta')
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('y in m (depth)','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
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');
%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(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(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
end
%
% exportfig(3, [fname parameter '.eps'],'bounds','tight', 'color','rgb', ...
% 'preview','none', 'resolution',200, 'lockaxes',1);
%3D plot of model
clear all;
nx=220; ny=100; nz=220; %ny:vertical
outx=1; outy=1; outz=1;
dh=0.15;
nx=nx/outx;ny=ny/outy;nz=nz/outz;
fignum=12;
%file_inp1='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/model/toy_real.vp';
%file_inp2='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/model/toy_real.vp';
file_inp1='../../par/model/rheinstetten.vs_it0';
file_inp2='../../par/model/rheinstetten.vs_it60';
phi1=0; %rotation angles to x-z plane of first and second plane
phi2=90;
rotaxes=[1,0,0]; %direction rotation axes [0,1,0] rotation of plane around vertical axis
% [1,0,0] rotation of plane around x-axis
%rotpoint=[120.8 140.8 80.8] ; %turning point [x y z] in meter
rotpoint=[(nx/4)*dh (15)*dh (3*nz/4)*dh] ;
%rotpoint=[56 60 60] ;
rotpoint2=[110 130 70] ;
%viewpoint=[122,26];
viewpoint=[-20,12];
caxis_value_1=0;%vs
caxis_value_2=300;%vs
%caxis_value_1=0;%vp
%caxis_value_2=2000;%vp
xcut1=1; xcut2=nx;
ycut1=1; ycut2=ny;
zcut1=1; zcut2=nz;
fid_rot=fopen(file_inp1,'r','ieee-le');
rot1=zeros(ny/outy,nx/outx,nz/outz);
rot1=fread(fid_rot,(nx*ny*nz),'float');
% fid_div=fopen(file_inp2,'r','ieee-le');
% div1=zeros(nz/outz,nx/outx,ny/outy);
% div1=fread(fid_div,(nx*ny*nz),'float');
%rot1=-rot1./max(max(rot1));
%rot1=rot1-div1;
%max(max(rot1))
%rot1=1./((rot1+3e-10));
%rot1=rot1.*div1*50000;
%rot1=-rot1./max(max(abs(rot1)));
%rot1=log10(rot1);
rot=reshape(rot1,ny/outy,nx/outx,nz/outz);
rot=rot(ycut1:ycut2,xcut1:xcut2,zcut1:zcut2);
nx=xcut2-xcut1+1;
ny=ycut2-ycut1+1;
nz=zcut2-zcut1+1;
%size(rot)
xp1=xcut1*dh; xp2=xcut2*dh; yp1=ycut1*dh; yp2=ycut2*dh; zp1=zcut1*dh; zp2=zcut2*dh;
x=xp1:dh*outx:xp2*outx;
y=yp1:dh*outy:yp2*outy;
z=zp1:dh*outz:zp2*outz;
%x=xcut1:dh*outx:xcut2*outx;
%y=ycut1:dh*outy:ycut2*outy;
%z=zcut1:dh*outz:zcut2*outz;
figure(12)
[Z,X,Y]=meshgrid(z,x,y);
xmin = min(X(:)); ymin = min(Y(:)); zmin = min(Z(:));
xmax = max(X(:)); ymax = max(Y(:)); zmax = max(Z(:));
%zplane=zeros(nx,nz);
%zplane(:,:)=rotpoint2(2);
%hslicez = surf(z,x,zplane);
%rotate(hslicez,rotaxes,phi1,rotpoint);
%xd1 = get(hslicez,'XData');
%yd1 = get(hslicez,'YData');
%zd1 = get(hslicez,'ZData');
%zplane=zeros(ny,nx);
%zplane(:,:)=rotpoint(3);
%hslicez = surf(x,y,zplane);
%rotate(hslicez,rotaxes,phi1,rotpoint);
%xd = get(hslicez,'XData');
%yd = get(hslicez,'YData');
%zd = get(hslicez,'ZData');
xd2 = rotpoint(1);
yd2 = rotpoint(2);
zd2 = rotpoint(3);
%zplane=zeros(nz,ny);
%zplane(:,:)=rotpoint(3);
%hslicez = surf(y,z,zplane);
%rotate(hslicez,rotaxes,phi2,rotpoint);
%xd2 = get(hslicez,'XData');
%yd2 = get(hslicez,'YData');
%zd2 = get(hslicez,'ZData');
%zplane=zeros(ny,nx);
%zplane(:,:)=rotpoint(3);
%hslicez = surf(x,y,zplane);
%rotate(hslicez,rotaxes,phi2,rotpoint);
%xd2 = get(hslicez,'XData');
%yd2 = get(hslicez,'YData');
%zd2 = get(hslicez,'ZData');
zplane=zeros(ny,nx);
zplane(:,:)=rotpoint(3);
hslicez = surf(x,y,zplane);
rotate(hslicez,[0 1 0],phi2,rotpoint);
xd3 = get(hslicez,'XData');
yd3 = get(hslicez,'YData');
zd3 = get(hslicez,'ZData');
%zplane=zeros(ny,nx);
%zplane(:,:)=rotpoint2(3);
%hslicez = surf(x,y,zplane);
%rotate(hslicez,[1 0 0],90,rotpoint2);
%xd4 = get(hslicez,'XData');
%yd4 = get(hslicez,'YData');
%zd4 = get(hslicez,'ZData');
rot=permute(rot,[2,3,1]);
%rott=size(rot)
%size(rot)
%size(Z)
figure(fignum)
%h3 = slice(Z,X,Y,rot,zd1,xd1,yd1);set(h3,'FaceColor','interp',...
% 'EdgeColor','none',...
% 'DiffuseStrength',.8)
% alpha(.5)
%hold on
%h = slice(Z,X,Y,rot,zd,xd,yd);set(h,'FaceColor','interp',...
% 'EdgeColor','none',...
% 'DiffuseStrength',.8)
%hold on
h2 = slice(Z,X,Y,rot,zd2,xd2,yd2);set(h2,'FaceColor','interp',...
'EdgeColor','none',...
'DiffuseStrength',.8)
hold on
% h2 = slice(Z,X,Y,rot,zd3,xd3,yd3);set(h2,'FaceColor','interp',...
% 'EdgeColor','none',...
% 'DiffuseStrength',.8)
hold on
plot3([rotpoint(3) rotpoint(3)],[xcut1*dh xcut1*dh],[ycut1*dh ycut2*dh],'-black','LineWidth',0.5);
plot3([rotpoint(3) rotpoint(3)],[xcut1*dh xcut2*dh],[ycut2*dh ycut2*dh],'-black','LineWidth',0.5);
plot3([rotpoint(3) rotpoint(3)],[xcut2*dh xcut2*dh],[ycut1*dh ycut2*dh],'-black','LineWidth',0.5);
plot3([rotpoint(3) rotpoint(3)],[xcut1*dh xcut2*dh],[ycut1*dh ycut1*dh],'-black','LineWidth',0.5);
hold on
plot3([zcut1*dh zcut2*dh],[xcut1*dh xcut1*dh],[rotpoint(2) rotpoint(2)],'-black','LineWidth',0.5);
plot3([zcut2*dh zcut2*dh],[xcut1*dh xcut2*dh],[rotpoint(2) rotpoint(2)],'-black','LineWidth',0.5);
plot3([zcut1*dh zcut2*dh],[xcut2*dh xcut2*dh],[rotpoint(2) rotpoint(2)],'-black','LineWidth',0.5);
plot3([zcut1*dh zcut1*dh],[xcut1*dh xcut2*dh],[rotpoint(2) rotpoint(2)],'-black','LineWidth',0.5);
hold on
plot3([zcut1*dh zcut2*dh],[rotpoint(1) rotpoint(1)],[ycut1*dh ycut1*dh],'-black','LineWidth',0.5);
plot3([zcut1*dh zcut2*dh],[rotpoint(1) rotpoint(1)],[ycut2*dh ycut2*dh],'-black','LineWidth',0.5);
plot3([zcut1*dh zcut1*dh],[rotpoint(1) rotpoint(1)],[ycut1*dh ycut2*dh],'-black','LineWidth',0.5);
plot3([zcut2*dh zcut2*dh],[rotpoint(1) rotpoint(1)],[ycut1*dh ycut2*dh],'-black','LineWidth',0.5);
hold on
plot3([zcut1*dh zcut2*dh],[rotpoint(1) rotpoint(1)],[rotpoint(2) rotpoint(2)],'k--','LineWidth',0.5);
plot3([rotpoint(3) rotpoint(3)],[xcut1*dh xcut2*dh],[rotpoint(2) rotpoint(2)],'k--','LineWidth',0.5);
plot3([rotpoint(3) rotpoint(3)],[rotpoint(1) rotpoint(1)],[ycut1*dh ycut2*dh],'k--','LineWidth',0.5);
% hold on
% hold on
% for(i=1:13)
% for(kk=1:13)
% plot3([20*dh+(kk-1)*10*dh 20*dh+(kk-1)*10*dh],[20*dh+(i-1)*10*dh 20*dh+(i-1)*10*dh],[24 24],'+blue','MarkerSize',3,'LineWidth',2.0);
% end
% end
% hold on
hold off
xlabel('z in m');
ylabel('x in m');
zlabel('y in m');
set(gca, 'xDir','normal')
set(gca, 'yDir','normal')
set(gca, 'zDir','reverse')
set(get(gca,'Ylabel'),'FontSize',14);
set(get(gca,'Ylabel'),'FontWeight','normal');
set(get(gca,'Xlabel'),'FontSize',14);
set(get(gca,'Xlabel'),'FontWeight','normal');
set(get(gca,'Zlabel'),'FontSize',14);
set(get(gca,'Zlabel'),'FontWeight','normal');
set(gca,'FontSize',14);
set(gca,'FontWeight','normal');
set(gca,'Linewidth',1.0);
colormap('jet')
cb = colorbar('vert');
xlabel(cb, 'v_p in m/s');
%set(cb,'YTick',[5800:200:6600]); % v_p
set(gca,'FontSize',14);
%load('MyColormapsgrad','mycmap');
% set(figure(fignum),'Colormap',mycmap)
view(viewpoint);
xlim([xcut1*dh xcut2*dh]);
ylim([zcut1*dh zcut2*dh]);
zlim([ycut1*dh ycut2*dh]);
caxis([caxis_value_1 caxis_value_2])
daspect([1,1,1]);
% axis tight
box on
% exportfig(fignum, 'toy_model_real_vp.eps','bounds','tight', 'color','rgb', ...
% 'preview','none', 'resolution',200, 'lockaxes',1);
......@@ -205,7 +205,7 @@ caxis([caxis_value_1 caxis_value_2])
set(gca,'ydir','normal');
xlabel('x in m','FontSize',fontsize)
ylabel('z in m','FontSize',fontsize)
title(['true ' parameter '-model. Slize at x=' num2str(sly) '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
......@@ -236,7 +236,7 @@ caxis([caxis_value_1 caxis_value_2])
set(gca,'ydir','normal');
xlabel('x in m','FontSize',fontsize)
ylabel('z in m','FontSize',fontsize)
title([parameter '-model. Iteration ' num2str(iteration) ' Slize at x=' num2str(sly) '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
......@@ -267,7 +267,7 @@ set(gca,'ydir','normal');
%set(gca,'xdir','reverse');
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['true ' parameter '-model. Slize at x=' num2str(slz) '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
......@@ -298,7 +298,7 @@ set(gca,'ydir','normal');
%set(gca,'xdir','reverse');
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title([parameter '-model. Iteration ' num2str(iteration) ' Slize at x=' num2str(slz) '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
......
#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
#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
#3
0.0000 1.7030 0.3569 0.0917
9.0000 2.0010 1.7450 0.2751
9.0000 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.5",
"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"
}
{
"MODELING PARAMETERS" : "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.5",
"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",
<