Commit 608c9d4a authored by laura.gassner's avatar laura.gassner

Merge branch 'development' of git.scc.kit.edu:GPIAG-Software/DENISE into development

parents ac564950 1b501c5d
close all
clear all
% This m-file is for making movies of wave propagation.
%%%%%%%%%%%%%%%%%%%%%%%% INPUT PARAMETER %%%%%%%%%%%%%%%%%%%%%%%%%
% Here the basic name of the binary snapshot files must be given:
% (The default extension is *.bin. The increasing number of the snapshot
% is added to the basic name, e.g. if ten snapshots were computed
% and the basic name is snap, then the filenames for the snapshots
% are snap1.bin, snap2.bin, ..., snap10.bin.
%filediv='snap/Khang_hom.bin.p.00';
% Receiver Positions
%yrec=3.6:0.66:69.6;
%xrec=45.5.*ones(length(yrec),1);
clear all
close all
% plot distance of sources and receivers
dnrec=10;
dnsour=1;
jseafloor=440;
% write model output files
WRITEMODE=1;
READ=1;
xs = 10.0;
ys = 0.48;
%rec=load('receiver_spheres.dat');
%xrec1=rec(:,1);
%yrec1=rec(:,2);
%xrec = xrec(1:dnrec:length(xrec));
%yrec = yrec(1:dnrec:length(yrec));
%source=load('sources_plot.dat')
%xshot=source(:,1);
%yshot=source(:,3);
%xshot = xshot(1:dnsour:length(xshot));
%yshot = yshot(1:dnsour:length(yshot));
% -------------------------------------------------------------------------
% P-Wave Velocity
% -------------------------------------------------------------------------
SMOOTH=0;
IMP=1;
vpmod = [1500.0 2000.0 3000.0];
vsmod = [866.0 1155.0 1732.0];
rhomod = [1929.0 2073.0 2294.0]
IDX=1; IDY=1;
IDXI=10; IDYI=10;
dh=0.1;
% gridsize and grid spacing (as specified in parameter-file)
NX1=1; NX2=500;
NY1=1; NY2=150;
%caxis_value1 = -1e-1;
%caxis_value2 = 1e-1;
caxis_value_vp1 = -1e-15;
caxis_value_vp2 = 1e-15;
%caxis_value_vp1 = -1e-2;
%caxis_value_vp2 = 1e-2;
caxis_value_vs1 = 0;
caxis_value_vs2 = 10;
caxis_value_rho1 = -1e-1;
caxis_value_rho2 = 1e-1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% grid size
nx=NX2-NX1+1
ny=NY2-NY1+1
if(READ==1)
clear vp;
clear vs;
clear rho;
% load model
file='jacobian_Test_p.old';
disp([' loading file ' file]);
fid=fopen(file,'r','ieee-le');
grad=fread(fid,[ny,nx],'float');
fclose(fid);
% load model
file='jacobian_Test_hessian';
disp([' loading file ' file]);
fid=fopen(file,'r','ieee-le');
hess=fread(fid,[ny,nx],'float');
fclose(fid);
end
%vp1 = vp(1:10:ny,1:10:nx);
%clear vp;
%vp = vp1;
%clear vp1;
%vs1 = vs(1:10:ny,1:10:nx);
%clear vs;
%vs = vs1;
%clear vs1;
%rho1 = rho(1:10:ny,1:10:nx);
%clear rho;
%rho = rho1;
%clear rho1;
NX1=1; NX2=500;
NY1=1; NY2=150;
dh=20.0;
% grid size
nx=NX2-NX1+1
ny=NY2-NY1+1
% plot range and increment
xp1=NX1*dh; xp2=NX2*dh; yp1=NY1*dh; yp2=NY2*dh;
% Computing range for axis and subscript range for the movie
x=xp1:IDX*dh:xp2;
y=yp1:IDY*dh:yp2;
epsilon = 5;
maxhess = max(max(hess));
for i=1:nx
for j=1:ny
depth = j.*dh;
ihess(j,i) = depth/(hess(j,i)+epsilon);
end
end
% apply damping of H^-1 in the PML layers
DAMPING=99;
amp=1.0-DAMPING./100.0; % amplitude at the edge of the numerical grid
ifw=10; % frame width in gridpoints
a=sqrt(-log(amp)./(ifw.^2));
for i=1:ifw
coeff(i)=exp(-(a.^2.*(ifw-i).^2));
%coeff(i) = 0.0;
end
% initialize array of coefficients with one
for j=1:ny
for i=1:nx
absorb_coeff(j,i)=1.0;
end
end
% coefficients for left and right boundary
yb=1;
ye=ny;
for i=1:ifw
ye=ny-i+1;
for j=yb:ye
absorb_coeff(j,i)=coeff(i);
end
end
yb=1;
ye=ny;
for i=1:ifw
ii=nx-i+1;
ye=ny-i+1;
for j=yb:ye
absorb_coeff(j,ii)=coeff(i);
end
end
% coefficients for bottom boundary
xb=1;
xe=nx;
for j=1:ifw
xb=j;
jj=ny-j+1;
xe=nx-j+1;
for i=xb:xe
absorb_coeff(jj,i)=coeff(j);
end
end
% apply damping in PML boundaries
ihess = ihess .* absorb_coeff;
subplot 311
imagesc(x,y,grad);
hold on;
plot(xs,ys,'k+');
%caxis([0 1]);
%caxis([caxis_value_vp1 caxis_value_vp2]);
%colorbar;
%set(gca,'YDir','normal');
%colorbar;
%axis equal;
%set(gca,'DataAspectRatio',[1 1 1]);
set(get(gca,'title'),'FontSize',12);
set(get(gca,'title'),'FontWeight','bold');
set(get(gca,'Ylabel'),'FontSize',12);
set(get(gca,'Ylabel'),'FontWeight','bold');
set(get(gca,'Xlabel'),'FontSize',12);
set(get(gca,'Xlabel'),'FontWeight','bold');
set(gca,'FontSize',12);
set(gca,'FontWeight','bold');
set(gca,'Box','on');
set(gca,'Linewidth',1.0);
axis ij
axis equal
axis tight
ylabel('y [m]');
iter_text=['Gradient'];
title(iter_text);
subplot 312
imagesc(x,y,ihess);
hold on;
plot(xs,ys,'k+');
%caxis([caxis_value_vs1 caxis_value_vs2]);
%colorbar;
%set(gca,'YDir','normal');
%colorbar;
%axis equal;
%set(gca,'DataAspectRatio',[1 1 1]);
set(get(gca,'title'),'FontSize',12);
set(get(gca,'title'),'FontWeight','bold');
set(get(gca,'Ylabel'),'FontSize',12);
set(get(gca,'Ylabel'),'FontWeight','bold');
set(get(gca,'Xlabel'),'FontSize',12);
set(get(gca,'Xlabel'),'FontWeight','bold');
set(gca,'FontSize',12);
set(gca,'FontWeight','bold');
set(gca,'Box','on');
set(gca,'Linewidth',1.0);
axis ij
axis equal
axis tight
ylabel('y [m]');
iter_text=['P (Hessian + 5I)^{-1}'];
title(iter_text);
%if(WRITEMODE==1)
%file1=['sphere_complex.vs'];
%fid1=fopen(file1,'w','ieee-le');
%fwrite(fid1,vs,'float')
%fclose(fid1);
%end
subplot 313
imagesc(x,y,ihess.*grad);
hold on;
plot(xs,ys,'k+');
% caxis([caxis_value_rho1 caxis_value_rho2]);
%colorbar;
%set(gca,'YDir','normal');
%colorbar;
%axis equal;
%set(gca,'DataAspectRatio',[1 1 1]);
set(get(gca,'title'),'FontSize',12);
set(get(gca,'title'),'FontWeight','bold');
set(get(gca,'Ylabel'),'FontSize',12);
set(get(gca,'Ylabel'),'FontWeight','bold');
set(get(gca,'Xlabel'),'FontSize',12);
set(get(gca,'Xlabel'),'FontWeight','bold');
set(gca,'FontSize',12);
set(gca,'FontWeight','bold');
set(gca,'Box','on');
set(gca,'Linewidth',1.0);
axis equal
axis tight
xlabel('x [m]');
ylabel('y [m]');
iter_text=['P (H + 5 I)^{-1} * Gradient'];
title(iter_text);
if(WRITEMODE==1)
file1=['taper.bin'];
fid1=fopen(file1,'w','ieee-le');
fwrite(fid1,ihess,'float')
fclose(fid1);
end
figure;
imagesc(x,y,ihess);
hold on;
plot(xs,ys,'k+');
% caxis([caxis_value_rho1 caxis_value_rho2]);
colorbar;
%set(gca,'YDir','normal');
%colorbar;
%axis equal;
%set(gca,'DataAspectRatio',[1 1 1]);
set(get(gca,'title'),'FontSize',12);
set(get(gca,'title'),'FontWeight','bold');
set(get(gca,'Ylabel'),'FontSize',12);
set(get(gca,'Ylabel'),'FontWeight','bold');
set(get(gca,'Xlabel'),'FontSize',12);
set(get(gca,'Xlabel'),'FontWeight','bold');
set(gca,'FontSize',12);
set(gca,'FontWeight','bold');
set(gca,'Box','on');
set(gca,'Linewidth',1.0);
axis equal
axis tight
xlabel('x [m]');
ylabel('y [m]');
iter_text=['H^{-1}'];
title(iter_text);
figure;
imagesc(x,y,absorb_coeff);
hold on;
plot(xs,ys,'k+');
% caxis([caxis_value_rho1 caxis_value_rho2]);
colorbar;
%set(gca,'YDir','normal');
%colorbar;
%axis equal;
%set(gca,'DataAspectRatio',[1 1 1]);
set(get(gca,'title'),'FontSize',12);
set(get(gca,'title'),'FontWeight','bold');
set(get(gca,'Ylabel'),'FontSize',12);
set(get(gca,'Ylabel'),'FontWeight','bold');
set(get(gca,'Xlabel'),'FontSize',12);
set(get(gca,'Xlabel'),'FontWeight','bold');
set(gca,'FontSize',12);
set(gca,'FontWeight','bold');
set(gca,'Box','on');
set(gca,'Linewidth',1.0);
axis equal
axis tight
xlabel('x [m]');
ylabel('y [m]');
iter_text=['damping coeff'];
title(iter_text);
%if(WRITEMODE==1)
%file1=['sphere_complex.rho'];
%fid1=fopen(file1,'w','ieee-le');
%fwrite(fid1,rho,'float')
%fclose(fid1);
%end
...@@ -28,6 +28,7 @@ endif ...@@ -28,6 +28,7 @@ endif
.PHONY: clean .PHONY: clean
clean: clean:
rm -rf ../bin/denise
$(MAKE) -C ../contrib/aff clean $(MAKE) -C ../contrib/aff clean
$(MAKE) -C ../contrib/fourier clean $(MAKE) -C ../contrib/fourier clean
$(MAKE) -C ../contrib/stfinv clean $(MAKE) -C ../contrib/stfinv clean
......
...@@ -44,7 +44,7 @@ ...@@ -44,7 +44,7 @@
"ACOUSTIC" : "0", "ACOUSTIC" : "0",
"Wavetype Computation" : "comment", "Wavetype Computation" : "comment",
"WAVETYPE" : "3", "WAVETYPE" : "2",
"Model" : "comment", "Model" : "comment",
"READMOD" : "0", "READMOD" : "0",
......
...@@ -5,25 +5,23 @@ ...@@ -5,25 +5,23 @@
############################################################### ###############################################################
# compiling all libraries and DENISE # compiling all libraries and DENISE
# make clean make clean
make denise MODEL=../genmod/toy_example_true.c make denise MODEL=../genmod/toy_example_true.c
pause
# starting DENISE for forward modeling # starting DENISE for forward modeling
# (depending on the MPI package you maybe have to adjust the following programme call, # (depending on the MPI package you maybe have to adjust the following programme call,
# e.g. when using openMPI you do not have to use the 'lamboot' command) # e.g. when using openMPI you do not have to use the 'lamboot' command)
#lamboot #lamboot
mpirun -np 4 ../bin/denise in_and_out/toy_example/toy_example_FW.json | tee in_and_out/toy_example/toy_example_FW.out mpirun -np 4 ../bin/denise in_and_out/toy_example/toy_example_FW_SH.json | tee in_and_out/toy_example/toy_example_FW_SH.out
# the forward modeled data have to be renamed for the inversion # the forward modeled data have to be renamed for the inversion
for (( i=1; i <= 5; i++ )) ; do for (( i=1; i <= 5; i++ )) ; do
mv su/measured_data/toy_example_vx.su.shot${i}.it1 su/measured_data/toy_example_x.su.shot${i} mv su/measured_data/toy_example_vx.su.shot${i}.it1 su/measured_data/toy_example_x.su.shot${i}
mv su/measured_data/toy_example_vy.su.shot${i}.it1 su/measured_data/toy_example_y.su.shot${i} mv su/measured_data/toy_example_vy.su.shot${i}.it1 su/measured_data/toy_example_y.su.shot${i}
mv su/measured_data/toy_example_vz.su.shot${i}.it1 su/measured_data/toy_example_z.su.shot${i} mv su/measured_data/toy_example_vz.su.shot${i}.it1 su/measured_data/toy_example_z.su.shot${i}
done done
exit
############################################################### ###############################################################
# running the inversion # # running the inversion #
......
...@@ -99,7 +99,7 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float ...@@ -99,7 +99,7 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
/* Debugging */ /* Debugging */
if(!ACOUSTIC) { if(!ACOUSTIC) {
sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iteration,w,POS[1],POS[2]); sprintf(jac,"%s_y_LBFGS_vs_it%d.bin.%i.%i",JACOBIAN,iteration,POS[1],POS[2]);
FP_JAC=fopen(jac,"wb"); FP_JAC=fopen(jac,"wb");
} }
...@@ -142,10 +142,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float ...@@ -142,10 +142,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
if(!ACOUSTIC) { if(!ACOUSTIC) {
fclose(FP_JAC); fclose(FP_JAC);
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin",JACOBIAN,iteration,w); sprintf(jac,"%s_y_LBFGS_vs_it%d.bin",JACOBIAN,iteration);
if (MYID==0) mergemod(jac,3); if (MYID==0) mergemod(jac,3);
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iteration,w,POS[1],POS[2]); sprintf(jac,"%s_y_LBFGS_vs_it%d.bin.%i.%i",JACOBIAN,iteration,POS[1],POS[2]);
remove(jac); remove(jac);
} }
......
...@@ -62,17 +62,17 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float * ...@@ -62,17 +62,17 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float *
if(w==0) w=N_LBFGS; if(w==0) w=N_LBFGS;
if(!ACOUSTIC){ if(!ACOUSTIC){
sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]); sprintf(jac,"%s_s_LBFGS_vs_it%d.bin.%i.%i",JACOBIAN,iter+1,POS[1],POS[2]);
FP_JAC=fopen(jac,"wb"); FP_JAC=fopen(jac,"wb");
} }
if(LBFGS_NPAR>1){ if(LBFGS_NPAR>1){
sprintf(jac2,"%s_s_LBFGS_rho_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]); sprintf(jac2,"%s_s_LBFGS_rho_it%d.bin.%i.%i",JACOBIAN,iter+1,POS[1],POS[2]);
FP_JAC2=fopen(jac2,"wb"); FP_JAC2=fopen(jac2,"wb");
} }
if(LBFGS_NPAR>2){ if(LBFGS_NPAR>2){
sprintf(jac2,"%s_s_LBFGS_vp_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]); sprintf(jac2,"%s_s_LBFGS_vp_it%d.bin.%i.%i",JACOBIAN,iter+1,POS[1],POS[2]);
FP_JAC3=fopen(jac2,"wb"); FP_JAC3=fopen(jac2,"wb");
} }
...@@ -344,30 +344,30 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float * ...@@ -344,30 +344,30 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float *
if(!ACOUSTIC){ if(!ACOUSTIC){
fclose(FP_JAC); fclose(FP_JAC);
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin",JACOBIAN,iter+1,w); sprintf(jac,"%s_s_LBFGS_vs_it%d.bin",JACOBIAN,iter+1);
if (MYID==0) mergemod(jac,3); if (MYID==0) mergemod(jac,3);
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]); sprintf(jac,"%s_s_LBFGS_vs_it%d.bin.%i.%i",JACOBIAN,iter+1,POS[1],POS[2]);
remove(jac); remove(jac);
} }
if(LBFGS_NPAR>1){ if(LBFGS_NPAR>1){
fclose(FP_JAC2); fclose(FP_JAC2);
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_s_LBFGS_rho_it%d_w%d.bin",JACOBIAN,iter+1,w); sprintf(jac,"%s_s_LBFGS_rho_it%d.bin",JACOBIAN,iter+1);
if (MYID==0) mergemod(jac,3); if (MYID==0) mergemod(jac,3);
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_s_LBFGS_rho_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]); sprintf(jac,"%s_s_LBFGS_rho_it%d.bin.%i.%i",JACOBIAN,iter+1,POS[1],POS[2]);
remove(jac); remove(jac);
} }
if(LBFGS_NPAR>2){ if(LBFGS_NPAR>2){
fclose(FP_JAC3); fclose(FP_JAC3);
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_s_LBFGS_vp_it%d_w%d.bin",JACOBIAN,iter+1,w); sprintf(jac,"%s_s_LBFGS_vp_it%d.bin",JACOBIAN,iter+1);
if (MYID==0) mergemod(jac,3); if (MYID==0) mergemod(jac,3);
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_s_LBFGS_vp_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]); sprintf(jac,"%s_s_LBFGS_vp_it%d.bin.%i.%i",JACOBIAN,iter+1,POS[1],POS[2]);
remove(jac); remove(jac);
} }
} }
......
...@@ -2584,16 +2584,16 @@ int main(int argc, char **argv){ ...@@ -2584,16 +2584,16 @@ int main(int argc, char **argv){
if(WAVETYPE==1 || WAVETYPE==3) { if(WAVETYPE==1 || WAVETYPE==3) {
We[j][i]=We[j][i]/We_max; We[j][i]=We[j][i]/We_max;
waveconv_shot[j][i] = waveconv_shot[j][i]/(We[j][i]*C_vp); waveconv_shot[j][i] = waveconv_shot[j][i]/(We[j][i]);
if(!ACOUSTIC){ if(!ACOUSTIC){
waveconv_u_shot[j][i] = waveconv_u_shot[j][i]/(We[j][i]*C_vs); waveconv_u_shot[j][i] = waveconv_u_shot[j][i]/(We[j][i]);
} }
waveconv_rho_shot[j][i] = waveconv_rho_shot[j][i]/(We[j][i]*C_rho); waveconv_rho_shot[j][i] = waveconv_rho_shot[j][i]/(We[j][i]);
} }
if(WAVETYPE==2 || WAVETYPE==3) { if(WAVETYPE==2 || WAVETYPE==3) {
We_SH[j][i]=We_SH[j][i]/We_max_SH; We_SH[j][i]=We_SH[j][i]/We_max_SH;
waveconv_u_shot_z[j][i] = waveconv_u_shot_z[j][i]/(We_SH[j][i]*C_vs); waveconv_u_shot_z[j][i] = waveconv_u_shot_z[j][i]/(We_SH[j][i]);
waveconv_rho_shot_z[j][i] = waveconv_rho_shot_z[j][i]/(We_SH[j][i]*C_rho); waveconv_rho_shot_z[j][i] = waveconv_rho_shot_z[j][i]/(We_SH[j][i]);
} }
} }
} }
...@@ -2714,11 +2714,11 @@ int main(int argc, char **argv){ ...@@ -2714,11 +2714,11 @@ int main(int argc, char **argv){
for (j=1;j<=NY;j=j+IDY){ for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){ for (i=1;i<=NX;i=i+IDX){
We_sum[j][i]=We_sum[j][i]/We_sum_max1; We_sum[j][i]=We_sum[j][i]/We_sum_max1;
waveconv[j][i] = waveconv[j][i]*We_sum[j][i]/C_vp; waveconv[j][i] = waveconv[j][i]*We_sum[j][i];
if(!ACOUSTIC){ if(!ACOUSTIC){
waveconv_u[j][i] = waveconv_u[j][i]*We_sum[j][i]/C_vs; waveconv_u[j][i] = waveconv_u[j][i]*We_sum[j][i];
} }
waveconv_rho[j][i] = waveconv_rho[j][i]*We_sum[j][i]/C_rho; waveconv_rho[j][i] = waveconv_rho[j][i]*We_sum[j][i];
} }
} }
} }
...@@ -2728,8 +2728,8 @@ int main(int argc, char **argv){ ...@@ -2728,8 +2728,8 @@ int main(int argc, char **argv){
for (j=1;j<=NY;j=j+IDY){ for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){ for (i=1;i<=NX;i=i+IDX){
We_sum_SH[j][i]=We_sum_SH[j][i]/We_sum_max1; We_sum_SH[j][i]=We_sum_SH[j][i]/We_sum_max1;
waveconv_u_z[j][i] = waveconv_u_z[j][i]*We_sum_SH[j][i]/C_vs; waveconv_u_z[j][i] = waveconv_u_z[j][i]*We_sum_SH[j][i];
waveconv_rho_z[j][i] = waveconv_rho_z[j][i]*We_sum_SH[j][i]/C_rho; waveconv_rho_z[j][i] = waveconv_rho_z[j][i]*We_sum_SH[j][i];
} }
} }
} }
......
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