Commit 7d089700 authored by laura.gassner's avatar laura.gassner

small adaptions for acoustic inversion, ac. TE now runs with L-BFGS

parent e7fc6d80
% This is a matlab script to plot the results of the acoustic toy example of the DENISE code.
clear all
% close all
clc
iteration_no=100;
% general informations about the model
nx=150; % grid points in x direction
ny=150; % grid points in y direction
DH=10.0; % spatial sampling interval
X=(1:1:nx)*DH; % definition of x axis for plots
Y=(1:1:ny)*DH; % definition of y axis for plots
FW=30; % width of PML frame in gridpoints
sources=[350.0, 550.0, 750.0, 950.0, 1150.0]; % x coordinate of source positions
% general settings for plots information
fontsize=12;
x_line=[FW FW nx-FW nx-FW]*DH;
y_line=[1 ny-FW ny-FW 1]*DH;
% reading the true model
true_vp=read_binary_matrix(nx,ny,['../par/model/mod_toy_example_ac_true_vp_it0.bin']);
true_rho=read_binary_matrix(nx,ny,['../par/model/mod_toy_example_ac_true_rho_it0.bin']);
% reading initial model
start_vp=read_binary_matrix(nx,ny,['../par/model/toy_example/mod_toy_example_ac_vp_it0.bin']);
start_rho=read_binary_matrix(nx,ny,['../par/model/toy_example/mod_toy_example_ac_rho_it0.bin']);
% reading inverted vp model
inv_vp=read_binary_matrix(nx,ny,['../par/model/toy_example/mod_toy_example_ac_vp_it' int2str(iteration_no) '.bin']);
% finding minimum and maximum values for colorbar
max_vp=max([max(true_vp(:)) max(start_vp(:)) max(inv_vp(:))]);
min_vp=min([min(true_vp(:)) min(start_vp(:)) min(inv_vp(:))]);
max_rho=max([max(true_rho(:)) max(start_rho(:))]);
min_rho=min([min(true_rho(:)) min(start_rho(:))]);
% Plot true model
figure('units','normalized','outerposition',[0.8 0 0.2 1])
subplot(3,2,1)
imagesc(X,Y,true_vp)
hold on
line(x_line,y_line,'LineStyle','--','Color','k')
plot(sources,10.0*ones(1,5),'*r','MarkerSize',12)
hold off
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['true P-wave velocity model in m/s'],'FontSize',fontsize)
colb=colorbar;
coll=get(colb,'ylabel');
set(coll,'String','vp in m/s');
caxis([min_vp max_vp])
axis equal tight
subplot(3,2,2)
imagesc(X,Y,true_rho)
hold on
line(x_line,y_line,'LineStyle','--','Color','k')
plot(sources,10.0*ones(1,5),'*r','MarkerSize',12)
hold off
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['true density model in kg/m^3'],'FontSize',fontsize)
colb=colorbar;
coll=get(colb,'ylabel');
set(coll,'String','\rho in kg/m^3');
caxis([min_rho max_rho])
axis equal tight
% Plot initial model
subplot(3,2,3)
imagesc(X,Y,start_vp)
hold on
line(x_line,y_line,'LineStyle','--','Color','k')
plot(sources,10.0*ones(1,5),'*r','MarkerSize',12)
hold off
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['initial P-wave velocity model in m/s'],'FontSize',fontsize)
colb=colorbar;
coll=get(colb,'ylabel');
set(coll,'String','vp in m/s');
caxis([min_vp max_vp])
axis equal tight
subplot(3,2,4)
imagesc(X,Y,start_rho)
hold on
line(x_line,y_line,'LineStyle','--','Color','k')
plot(sources,10.0*ones(1,5),'*r','MarkerSize',12)
hold off
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['initial density model in kg/m^3'],'FontSize',fontsize)
colb=colorbar;
coll=get(colb,'ylabel');
set(coll,'String','\rho in kg/m^3');
caxis([min_rho max_rho])
axis equal tight
% Plot obtained model
subplot(3,2,5)
imagesc(X,Y,inv_vp)
hold on
line(x_line,y_line,'LineStyle','--','Color','k')
plot(sources,10.0*ones(1,5),'*r','MarkerSize',12)
hold off
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['inverted P-wave velocity model in m/s',10,'(iteration step ' int2str(iteration_no) ')'],'FontSize',fontsize)
colb=colorbar;
coll=get(colb,'ylabel');
set(coll,'String','vp in m/s');
caxis([min_vp max_vp])
axis equal tight
% Plot vertical velocity provfiles for S-wave velocity model
x1=500;
x2=750;
subplot(3,2,6)
plot(true_vp(:,round(x1/DH)),Y,'Color',[0.7 0.7 0.7],'LineWidth',6,'LineStyle','-');
hold on
plot(start_vp(:,round(x1/DH)),Y,'k','LineWidth',3,'LineStyle','--');
plot(inv_vp(:,round(x1/DH)),Y,'b','LineWidth',3);
plot(inv_vp(:,round(x2/DH)),Y,'r','LineWidth',3,'LineStyle','--');
line([min_vp-10 max_vp+10],[(ny-FW)*DH (ny-FW)*DH],'LineStyle','--','Color','k')
hold off
xlim([min_vp-10 max_vp+10])
set(gca,'ydir','reverse')
xlabel('velocity in m/s','fontsize',fontsize)
ylabel('depth in m','fontsize',fontsize)
title('Vertical velocity profiles','fontsize',fontsize)
legend('true','initial',['x=' num2str(x1) 'm'],['x=' num2str(x2) 'm']);
hold off
grid on
......@@ -83,38 +83,47 @@
"QUELLTYPB" : "4",
"MISFIT_LOG_FILE" : "LOG_toy_example_ac.dat",
"Inversion for density" : "comment",
"INV_RHO_ITER" : "100",
"INV_VP_ITER" : "1",
"Output of inverted models" : "comment",
"INV_MODELFILE" : "model/toy_example/mod_toy_example_ac",
"nfstart" : "1",
"nf" : "1",
"Output of gradients" : "comment",
"JACOBIAN" : "jacobian/toy_example/jac_toy_example_ac",
"nfstart_jac" : "1",
"nf_jac" : "1",
"Gradient-Method" : "comment",
"GRAD_METHOD" : "1",
"GRAD_METHOD" : "2",
"LBFGS_STEP_LENGTH" : "1",
"N_LBFGS" : "5",
"N_LBFGS" : "10",
"Wolfe Condition" : "comment",
"WOLFE_CONDITION" : "1",
"WOLFE_NUM_TEST" : "5",
"WOLFE_TRY_OLD_STEPLENGTH" : "1",
"WOLFE_C1_SL" : "0.0",
"Approx. Hessian" : "comment",
"EPRECOND" : "0",
"EPRECOND" : "3",
"EPSILON_WE" : "0.005",
"EPRECOND_ITER" : "0",
"EPRECOND_PER_SHOT" : "0",
"Workflow" : "comment",
"USE_WORKFLOW" : "0",
"FILE_WORKFLOW" : "workflow.txt",
"Inversion for density" : "comment",
"INV_RHO_ITER" : "100",
"INV_VP_ITER" : "1",
"Step length estimation" : "comment",
"EPS_SCALE" : "0.01",
"STEPMAX" : "4",
"SCALEFAC" : "5.0",
"TESTSHOT_START , TESTSHOT_END , TESTSHOT_INCR" : "1 , 5 , 2",
"Output of inverted models" : "comment",
"INV_MODELFILE" : "model/toy_example/mod_toy_example_ac",
"nfstart" : "1",
"nf" : "1",
"Output of gradients" : "comment",
"JACOBIAN" : "jacobian/toy_example/jac_toy_example_ac",
"nfstart_jac" : "1",
"nf_jac" : "1",
"Gradient calculation" : "comment",
"LNORM" : "2",
......@@ -122,12 +131,6 @@
"NORMALIZE" : "0",
"DTINV" : "1",
"Step length estimation" : "comment",
"EPS_SCALE" : "0.01",
"STEPMAX" : "4",
"SCALEFAC" : "5.0",
"TESTSHOT_START , TESTSHOT_END , TESTSHOT_INCR" : "1 , 5 , 2",
"Termination of the programmme" : "comment",
"PRO" : "0.01",
......@@ -144,6 +147,15 @@
"TAPER_FILE_NAME_U" : "taper_te_ac_u.bin",
"TAPER_FILE_NAME_RHO" : "taper_te_ac_rho.bin",
"Definition of inversion for source time function" : "comment",
"INV_STF" : "0",
"PARA" : "fdlsq:tshift=0.0",
"N_STF" : "10",
"N_STF_START" : "1",
"TAPER_STF" : "1",
"TRKILL_STF" : "0",
"TRKILL_FILE_STF" : "./trace_kill/trace_kill.dat",
"Upper and lower limits for model parameters" : "comment",
"VPUPPERLIM" : "5000",
"VPLOWERLIM" : "0",
......@@ -156,15 +168,6 @@
"MODEL_FILTER" : "0",
"FILT_SIZE" : "3",
"Definition of inversion for source time function" : "comment",
"INV_STF" : "0",
"PARA" : "fdlsq:tshift=0.0",
"N_STF" : "10",
"N_STF_START" : "1",
"TAPER_STF" : "1",
"TRKILL_STF" : "0",
"TRKILL_FILE_STF" : "./trace_kill/trace_kill.dat",
"Frequency filtering during inversion" : "comment",
"TIME_FILT" : "0",
"F_HP" : "1",
......@@ -180,6 +183,7 @@
"Time windowing and damping" : "comment",
"TIMEWIN" : "0",
"TW_IND" : "1",
"PICKS_FILE" : "./picked_times/PickedTimes",
"TWLENGTH_PLUS" : "4.0",
"TWLENGTH_MINUS" : "0.01",
......
......@@ -2585,7 +2585,9 @@ int main(int argc, char **argv){
if(WAVETYPE==1 || WAVETYPE==3) {
We[j][i]=We[j][i]/We_max;
waveconv_shot[j][i] = waveconv_shot[j][i]/(We[j][i]*C_vp);
waveconv_u_shot[j][i] = waveconv_u_shot[j][i]/(We[j][i]*C_vs);
if(!ACOUSTIC){
waveconv_u_shot[j][i] = waveconv_u_shot[j][i]/(We[j][i]*C_vs);
}
waveconv_rho_shot[j][i] = waveconv_rho_shot[j][i]/(We[j][i]*C_rho);
}
if(WAVETYPE==2 || WAVETYPE==3) {
......@@ -2713,7 +2715,9 @@ int main(int argc, char **argv){
for (i=1;i<=NX;i=i+IDX){
We_sum[j][i]=We_sum[j][i]/We_sum_max1;
waveconv[j][i] = waveconv[j][i]*We_sum[j][i]/C_vp;
waveconv_u[j][i] = waveconv_u[j][i]*We_sum[j][i]/C_vs;
if(!ACOUSTIC){
waveconv_u[j][i] = waveconv_u[j][i]*We_sum[j][i]/C_vs;
}
waveconv_rho[j][i] = waveconv_rho[j][i]*We_sum[j][i]/C_rho;
}
}
......@@ -4223,11 +4227,11 @@ int main(int argc, char **argv){
if(WOLFE_CONDITION){
free_matrix(waveconv_old,-nd+1,NY+nd,-nd+1,NX+nd);
free_matrix(waveconv_u_old,-nd+1,NY+nd,-nd+1,NX+nd);
if(!ACOUSTIC) free_matrix(waveconv_u_old,-nd+1,NY+nd,-nd+1,NX+nd);
free_matrix(waveconv_rho_old,-nd+1,NY+nd,-nd+1,NX+nd);
free_matrix(waveconv_up,-nd+1,NY+nd,-nd+1,NX+nd);
free_matrix(waveconv_u_up,-nd+1,NY+nd,-nd+1,NX+nd);
if(!ACOUSTIC) free_matrix(waveconv_u_up,-nd+1,NY+nd,-nd+1,NX+nd);
free_matrix(waveconv_rho_up,-nd+1,NY+nd,-nd+1,NX+nd);
}
......
......@@ -190,7 +190,38 @@ void read_par_json(FILE *fp, char *fileinp){
}
}
if (get_int_from_objectlist("ACOUSTIC",number_readobjects,&ACOUSTIC,varname_list, value_list)){
ACOUSTIC=0;
fprintf(fp,"Variable ACOUSTIC is set to default value %d.\n",ACOUSTIC);}
if (get_int_from_objectlist("WAVETYPE",number_readobjects,&WAVETYPE,varname_list, value_list)){
WAVETYPE=1;
fprintf(fp,"Variable WAVETYPE is set to default value %d.\n",WAVETYPE);
} else {
if (ACOUSTIC && WAVETYPE!=1) {
WAVETYPE=1;
fprintf(fp,"For acoustic modelling WAVETYPE is set to %d.\n",WAVETYPE);
}
if(WAVETYPE==3) {
if (get_int_from_objectlist("JOINT_INVERSION_PSV_SH_TYPE",number_readobjects,&JOINT_INVERSION_PSV_SH_TYPE,varname_list, value_list)){
JOINT_INVERSION_PSV_SH_TYPE=1;
fprintf(fp,"Variable JOINT_INVERSION_PSV_SH_TYPE is set to default value %d.\n",JOINT_INVERSION_PSV_SH_TYPE);
} else {
/* Check herer possible dependencies */
}
if (get_float_from_objectlist("JOINT_INVERSION_PSV_SH_ALPHA_VS",number_readobjects,&JOINT_INVERSION_PSV_SH_ALPHA_VS,varname_list, value_list)){
JOINT_INVERSION_PSV_SH_ALPHA_VS=0.5;
fprintf(fp,"Variable JOINT_INVERSION_PSV_SH_ALPHA_VS is set to default value %f.\n",JOINT_INVERSION_PSV_SH_ALPHA_VS);
} else {
/* Check herer possible dependencies */
}
if (get_float_from_objectlist("JOINT_INVERSION_PSV_SH_ALPHA_RHO",number_readobjects,&JOINT_INVERSION_PSV_SH_ALPHA_RHO,varname_list, value_list)){
JOINT_INVERSION_PSV_SH_ALPHA_RHO=0.5;
fprintf(fp,"Variable JOINT_INVERSION_PSV_SH_ALPHA_RHO is set to default value %f.\n",JOINT_INVERSION_PSV_SH_ALPHA_RHO);
} else {
/* Check herer possible dependencies */
}
}
}
if (get_int_from_objectlist("QUELLART",number_readobjects,&QUELLART,varname_list, value_list))
err("Variable QUELLART could not be retrieved from the json input file!");
......@@ -412,44 +443,11 @@ void read_par_json(FILE *fp, char *fileinp){
fprintf(fp,"Reference frequency for viscoelastic modeling is set to center frequency of the source wavelet.\n");}
}
if (get_int_from_objectlist("ACOUSTIC",number_readobjects,&ACOUSTIC,varname_list, value_list)){
ACOUSTIC=0;
fprintf(fp,"Variable ACOUSTIC is set to default value %d.\n",ACOUSTIC);}
/*=================================
section inversion parameters
=================================*/
if (get_int_from_objectlist("WAVETYPE",number_readobjects,&WAVETYPE,varname_list, value_list)){
WAVETYPE=1;
fprintf(fp,"Variable WAVETYPE is set to default value %d.\n",WAVETYPE);
} else {
if (ACOUSTIC && WAVETYPE!=1) {
WAVETYPE=1;
fprintf(fp,"For acoustic modelling WAVETYPE is set to %d.\n",WAVETYPE);
}
if(WAVETYPE==3) {
if (get_int_from_objectlist("JOINT_INVERSION_PSV_SH_TYPE",number_readobjects,&JOINT_INVERSION_PSV_SH_TYPE,varname_list, value_list)){
JOINT_INVERSION_PSV_SH_TYPE=1;
fprintf(fp,"Variable JOINT_INVERSION_PSV_SH_TYPE is set to default value %d.\n",JOINT_INVERSION_PSV_SH_TYPE);
} else {
/* Check herer possible dependencies */
}
if (get_float_from_objectlist("JOINT_INVERSION_PSV_SH_ALPHA_VS",number_readobjects,&JOINT_INVERSION_PSV_SH_ALPHA_VS,varname_list, value_list)){
JOINT_INVERSION_PSV_SH_ALPHA_VS=0.5;
fprintf(fp,"Variable JOINT_INVERSION_PSV_SH_ALPHA_VS is set to default value %f.\n",JOINT_INVERSION_PSV_SH_ALPHA_VS);
} else {
/* Check herer possible dependencies */
}
if (get_float_from_objectlist("JOINT_INVERSION_PSV_SH_ALPHA_RHO",number_readobjects,&JOINT_INVERSION_PSV_SH_ALPHA_RHO,varname_list, value_list)){
JOINT_INVERSION_PSV_SH_ALPHA_RHO=0.5;
fprintf(fp,"Variable JOINT_INVERSION_PSV_SH_ALPHA_RHO is set to default value %f.\n",JOINT_INVERSION_PSV_SH_ALPHA_RHO);
} else {
/* Check herer possible dependencies */
}
}
}
if (get_int_from_objectlist("INVMAT1",number_readobjects,&INVMAT1,varname_list, value_list))
err("Variable INVMAT1 could not be retrieved from the json input file!");
else
......
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