Commit 14beb822 authored by fabian.kuehn's avatar fabian.kuehn

Added Reverse Time Migration (RTM)

parent 04ca7c45
......@@ -17,7 +17,7 @@
-----------------------------------------------------------------------------------------*/
/* ----------------------------------------------------------------------
* This is program IFOS Version 2.0.3
* This is program IFOS Version 2.0.1
* Inversion of Full Observerd Seismograms
*
* ----------------------------------------------------------------------*/
......@@ -53,14 +53,15 @@ int main(int argc, char **argv){
int sum_killed_traces=0, sum_killed_traces_testshots=0, killed_traces=0, killed_traces_testshots=0;
int *ptr_killed_traces=&killed_traces, *ptr_killed_traces_testshots=&killed_traces_testshots;
float energy, energy_sum, energy_all_shots, energy_sum_all_shots = 0.0;
float energy, energy_sum, energy_all_shots, energy_sum_all_shots;
float energy_SH, energy_sum_SH, energy_all_shots_SH, energy_sum_all_shots_SH;
float L2_SH, L2sum_SH, L2_all_shots_SH, L2sum_all_shots_SH;
// Pointer for dynamic wavefields:
float ** psxx, ** psxy, ** psyy, ** psxz, ** psyz, **psp, ** ux, ** uy, ** uxy, ** uyx, ** u, ** Vp0, ** uttx, ** utty, ** Vs0, ** Rho0;
float ** psxx, ** psxy, ** psyy, ** psxz, ** psyz, **psp, ** ux, ** uy, ** uxy, ** uyx, ** Vp0, ** uttx, ** utty, ** Vs0, ** Rho0;
float ** pvx, ** pvy, ** pvz, **waveconv, **waveconv_lam, **waveconv_mu, **waveconv_rho, **waveconv_rho_s, **waveconv_u, **waveconvtmp, **wcpart, **wavejac,**waveconv_rho_s_z,**waveconv_u_z,**waveconv_rho_z;
float **waveconv_shot, **waveconv_u_shot, **waveconv_rho_shot, **waveconv_u_shot_z, **waveconv_rho_shot_z;
float **waveconv_rtm_shot, **waveconv_rtm;
float ** pvxp1, ** pvyp1, ** pvzp1, ** pvxm1, ** pvym1, ** pvzm1;
float ** gradg, ** gradp,** gradg_rho, ** gradp_rho, ** gradg_u, ** gradp_u, ** gradp_u_z,** gradp_rho_z;
float ** prho,** prhonp1, **prip=NULL, **prjp=NULL, **pripnp1=NULL, **prjpnp1=NULL, ** ppi, ** pu, ** punp1, ** puipjp, ** ppinp1;
......@@ -157,15 +158,9 @@ int main(int argc, char **argv){
int nfrq=0;
int FREQ_NR=1;
float JOINT_EQUAL_PSV=0.0, JOINT_EQUAL_SH=0.0;
float JOINT_EQUAL_PSV_all=0.0, JOINT_EQUAL_SH_all=0.0;
int JOINT_EQUAL_new_max=1;
FILE *fprec, *FPL2;
FILE *FPL2_JOINT;
char L2_joint_log[STRING_SIZE];
/* General parameters */
int nt_out;
......@@ -244,13 +239,9 @@ int main(int argc, char **argv){
NX = IENDX;
NY = IENDY;
/* Reading source positions from SOURCE_FILE */
srcpos=sources(&nsrc);
nsrc_glob=nsrc;
ishot=0;
if (SEISMO){
recpos=receiver(&ntr, srcpos, ishot);
recpos=receiver(FP, &ntr);
recswitch = ivector(1,ntr);
recpos_loc = splitrec(recpos,&ntr_loc, ntr, recswitch);
ntr_glob=ntr;
......@@ -418,8 +409,6 @@ int main(int argc, char **argv){
uxz = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
uyz = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
}
}else{
u = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
}
switch (WAVETYPE) {
......@@ -528,6 +517,11 @@ int main(int argc, char **argv){
waveconv_lam = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_shot = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
if(RTM==1){
waveconv_rtm = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
waveconv_rtm_shot = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
}
waveconvtmp = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
wcpart = matrix(1,3,1,3);
wavejac = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
......@@ -715,13 +709,104 @@ int main(int argc, char **argv){
}
if (ntr>0){
alloc_sections(ntr,ns,&sectionvx,&sectionvy,&sectionvz,&sectionp,&sectionpnp1,&sectionpn,&sectioncurl,&sectiondiv,
&sectionpdata,&sectionpdiff,&sectionpdiffold,&sectionvxdata,&sectionvxdiff,&sectionvxdiffold,&sectionvydata,
&sectionvydiff,&sectionvydiffold,&sectionvzdata,&sectionvzdiff,&sectionvzdiffold);
switch (SEISMO){
case 1 : /* particle velocities only */
switch (WAVETYPE) {
case 1:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
sectionvz=matrix(1,ntr,1,ns);
break;
}
break;
case 2 : /* pressure only */
sectionp=matrix(1,ntr,1,ns);
sectionpnp1=matrix(1,ntr,1,ns);
sectionpn=matrix(1,ntr,1,ns);
break;
case 3 : /* curl and div only */
sectioncurl=matrix(1,ntr,1,ns);
sectiondiv=matrix(1,ntr,1,ns);
break;
case 4 : /* everything */
switch (WAVETYPE) {
case 1:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
sectionvz=matrix(1,ntr,1,ns);
break;
}
sectioncurl=matrix(1,ntr,1,ns);
sectiondiv=matrix(1,ntr,1,ns);
sectionp=matrix(1,ntr,1,ns);
break;
case 5 : /* everything except curl and div*/
switch (WAVETYPE) {
case 1:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
sectionvz=matrix(1,ntr,1,ns);
break;
}
sectionp=matrix(1,ntr,1,ns);
break;
}
}
/* Memory for seismic data */
sectionread=matrix(1,ntr_glob,1,ns);
sectionpdata=matrix(1,ntr,1,ns);
sectionpdiff=matrix(1,ntr,1,ns);
sectionpdiffold=matrix(1,ntr,1,ns);
switch (WAVETYPE) {
case 1:
sectionvxdata=matrix(1,ntr,1,ns);
sectionvxdiff=matrix(1,ntr,1,ns);
sectionvxdiffold=matrix(1,ntr,1,ns);
sectionvydata=matrix(1,ntr,1,ns);
sectionvydiff=matrix(1,ntr,1,ns);
sectionvydiffold=matrix(1,ntr,1,ns);
break;
case 2:
sectionvzdata=matrix(1,ntr,1,ns);
sectionvzdiff=matrix(1,ntr,1,ns);
sectionvzdiffold=matrix(1,ntr,1,ns);
break;
case 3:
sectionvxdata=matrix(1,ntr,1,ns);
sectionvxdiff=matrix(1,ntr,1,ns);
sectionvxdiffold=matrix(1,ntr,1,ns);
sectionvydata=matrix(1,ntr,1,ns);
sectionvydiff=matrix(1,ntr,1,ns);
sectionvydiffold=matrix(1,ntr,1,ns);
sectionvzdata=matrix(1,ntr,1,ns);
sectionvzdiff=matrix(1,ntr,1,ns);
sectionvzdiffold=matrix(1,ntr,1,ns);
break;
}
/* Memory for inversion for source time function */
if((INV_STF==1)||(TIME_FILT==1) || (TIME_FILT==2)){
......@@ -770,6 +855,10 @@ int main(int argc, char **argv){
MPI_Barrier(MPI_COMM_WORLD);
/* Reading source positions from SOURCE_FILE */
srcpos=sources(&nsrc);
nsrc_glob=nsrc;
if(FORWARD_ONLY==0&&USE_WORKFLOW){
read_workflow(FILE_WORKFLOW,&workflow, &workflow_lines,workflow_header);
}
......@@ -867,8 +956,6 @@ int main(int argc, char **argv){
gradient_optimization=1;
}
/* Reset fail status of parabolic step length search */
step3=0;
}
if (MYID==0){
......@@ -957,8 +1044,8 @@ int main(int argc, char **argv){
if(PARAMETERIZATION==2){
Vp0[j][i] = sqrt((ppi[j][i]+2.0*pu[j][i])/prho[j][i]);
Vs0[j][i] = sqrt((pu[j][i])/prho[j][i]);
Vp0[j][i] = sqrt((ppi[j][i]+2.0*pu[j][i])*prho[j][i]);
Vs0[j][i] = sqrt((pu[j][i])*prho[j][i]);
Rho0[j][i] = prho[j][i];
}
......@@ -987,41 +1074,27 @@ int main(int argc, char **argv){
C_rho = rho_avg*rho_avg;
}
/* Seperate PSV and SH logging in case of a joint inversion */
if(WAVETYPE==3){
sprintf(L2_joint_log,"%s_JOINT",MISFIT_LOG_FILE);
}
/* Open Log File for L2 norm */
if(!FORWARD_ONLY && MYID==0){
if(FORWARD_ONLY!=1){
if(MYID==0){
if(iter==1){
FPL2=fopen(MISFIT_LOG_FILE,"w");
/* Write header for misfit log file */
if(GRAD_METHOD==1&&VERBOSE) {
if (TIME_FILT==0){
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t GAMMA \n");}
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \n");}
else{
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t F_LOW_PASS \t GAMMA \n");
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t F_LOW_PASS \n");
}
}
if(WAVETYPE==3) FPL2_JOINT=fopen(L2_joint_log,"w");
} else {
FPL2=fopen(MISFIT_LOG_FILE,"a");
if(WAVETYPE==3) FPL2_JOINT=fopen(L2_joint_log,"a");
}
if(iter>1){FPL2=fopen(MISFIT_LOG_FILE,"a");}
}
}
/* initialization of L2 calculation */
L2=0.0;
Lcount=0;
energy=0.0;
L2_all_shots=0.0;
energy_all_shots=0.0;
......@@ -1044,6 +1117,7 @@ int main(int argc, char **argv){
for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){
waveconv[j][i]=0.0;
if(RTM==1) waveconv_rtm[j][i]=0.0;
waveconv_rho[j][i]=0.0;
if(!ACOUSTIC) waveconv_u[j][i]=0.0;
}
......@@ -1086,28 +1160,6 @@ int main(int argc, char **argv){
for (ishot=1;ishot<=nshots;ishot+=SHOTINC){
if (SEISMO && READREC==2){
if (ntr>0) {
dealloc_sections(ntr,ns,recpos_loc,sectionvx,sectionvy,sectionvz,sectionp,sectionpnp1,sectionpn,sectioncurl,sectiondiv,
sectionpdata,sectionpdiff,sectionpdiffold,sectionvxdata,sectionvxdiff,sectionvxdiffold,sectionvydata,
sectionvydiff,sectionvydiffold,sectionvzdata,sectionvzdiff,sectionvzdiffold);
}
free_imatrix(recpos,1,3,1,ntr_glob);
recpos=receiver(&ntr, srcpos, ishot);
recpos_loc = splitrec(recpos,&ntr_loc, ntr, recswitch);
ntr_glob=ntr;
ntr=ntr_loc;
if (ntr>0){
alloc_sections(ntr,ns,&sectionvx,&sectionvy,&sectionvz,&sectionp,&sectionpnp1,&sectionpn,&sectioncurl,&sectiondiv,
&sectionpdata,&sectionpdiff,&sectionpdiffold,&sectionvxdata,&sectionvxdiff,&sectionvxdiffold,&sectionvydata,
&sectionvydiff,&sectionvydiffold,&sectionvzdata,&sectionvzdiff,&sectionvzdiffold);
}
if (ntr) group_id = 1;
else group_id = 0;
MPI_Comm_split(MPI_COMM_WORLD, group_id, MYID, &MPI_COMM_NTR);
MPI_Comm_rank(MPI_COMM_NTR, &myid_ntr);
}
SOURCE_SHAPE = SOURCE_SHAPE_OLD;
if(WAVETYPE==2 || WAVETYPE==3) SOURCE_SHAPE_SH=SOURCE_SHAPE_OLD_SH;
......@@ -1258,7 +1310,7 @@ int main(int argc, char **argv){
if(!ACOUSTIC) {
update_s_elastic_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, ppi, pu, puipjp, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
} else {
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, u, ppi, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, ppi, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
}
if (WAVETYPE==2 || WAVETYPE==3) {
......@@ -1339,18 +1391,19 @@ int main(int argc, char **argv){
if(LNORM==8){
calc_envelope(fulldata_vy,fulldata_vy,ns,ntr_glob);
calc_envelope(fulldata_vx,fulldata_vx,ns,ntr_glob);}
// if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
if (MYID==0){
saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);}
break;
case 2 : /* pressure only */
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD);
// if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
break;
case 3 : /* curl and div only */
catseis(sectiondiv, fulldata_div, recswitch, ntr_glob, MPI_COMM_WORLD);
catseis(sectioncurl, fulldata_curl, recswitch, ntr_glob, MPI_COMM_WORLD);
// if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
break;
case 4 : /* everything */
......@@ -1364,7 +1417,7 @@ int main(int argc, char **argv){
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD);
catseis(sectiondiv, fulldata_div, recswitch, ntr_glob, MPI_COMM_WORLD);
catseis(sectioncurl, fulldata_curl, recswitch, ntr_glob, MPI_COMM_WORLD);
// if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
break;
case 5 : /* everything except curl and div*/
......@@ -1376,7 +1429,7 @@ int main(int argc, char **argv){
catseis(sectionvz, fulldata_vz, recswitch, ntr_glob, MPI_COMM_WORLD);
}
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD);
// if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
break;
} /* end of switch (SEISMO) */
......@@ -1572,6 +1625,7 @@ int main(int argc, char **argv){
for(j=1;j<=NY;j=j+IDY){
for(i=1;i<=NX;i=i+IDX){
waveconv_shot[j][i]=0.0;
if(RTM==1) waveconv_rtm_shot[j][i]=0.0;
waveconv_rho_shot[j][i]=0.0;
}
}
......@@ -1700,7 +1754,7 @@ int main(int argc, char **argv){
if(!ACOUSTIC) {
update_s_elastic_PML(1, NX, 1, NY, pvx, pvy, ux, uy, uxy, uyx, psxx, psyy, psxy, ppi, pu, puipjp, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
} else {
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, u, ppi, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, ppi, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
}
if (WAVETYPE==2 || WAVETYPE==3) {
......@@ -1790,11 +1844,7 @@ int main(int argc, char **argv){
}else{
for (j=1;j<=NY;j=j+IDYI){
for (i=1;i<=NX;i=i+IDXI){
if(VELOCITY==0){
forward_prop_p[j][i][hin]=psp[j][i];
}else{
forward_prop_p[j][i][hin]=u[j][i];
}
}
}
}
......@@ -1870,7 +1920,8 @@ int main(int argc, char **argv){
if(LNORM==8){
calc_envelope(fulldata_vy,fulldata_vy,ns,ntr_glob);
calc_envelope(fulldata_vx,fulldata_vx,ns,ntr_glob);}
if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
if (MYID==0){
saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);}
break;
case 2 : /* pressure only */
......@@ -1951,8 +2002,8 @@ int main(int argc, char **argv){
}
L2=calc_res(sectionvxdata,sectionvx,sectionvxdiff,sectionvxdiffold,ntr,ns,LNORM,L2,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
if(swstestshot==1){energy=calc_energy(sectionvxdata,ntr,ns,energy, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);}
energy_all_shots=calc_energy(sectionvxdata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
L2_all_shots=calc_misfit(sectionvxdata,sectionvx,ntr,ns,LNORM,L2_all_shots,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
energy_all_shots=calc_energy(sectionvxdata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
/*fprintf(FP,"Energy vxdata for PE %d: %f\n\n", MYID,energy);*/
} /* end ADJOINT_TYPE */
......@@ -1974,8 +2025,8 @@ int main(int argc, char **argv){
}
L2=calc_res(sectionvydata,sectionvy,sectionvydiff,sectionvydiffold,ntr,ns,LNORM,L2,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
if(swstestshot==1){energy=calc_energy(sectionvydata,ntr,ns,energy, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);}
energy_all_shots=calc_energy(sectionvydata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
L2_all_shots=calc_misfit(sectionvydata,sectionvy,ntr,ns,LNORM,L2_all_shots,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
energy_all_shots=calc_energy(sectionvydata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
/*fprintf(FP,"Energy vydata for PE %d: %f\n\n", MYID,energy); */
} /* end ADJOINT_TYPE */
......@@ -1996,8 +2047,8 @@ int main(int argc, char **argv){
}
L2=calc_res(sectionpdata,sectionp,sectionpdiff,sectionpdiffold,ntr,ns,LNORM,L2,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
if(swstestshot==1){energy=calc_energy(sectionpdata,ntr,ns,energy, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);}
energy_all_shots=calc_energy(sectionpdata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
L2_all_shots=calc_misfit(sectionpdata,sectionp,ntr,ns,LNORM,L2_all_shots,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
energy_all_shots=calc_energy(sectionpdata,ntr,ns,energy_all_shots, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
} /* end ADJOINT_TYPE */
}
......@@ -2018,8 +2069,8 @@ int main(int argc, char **argv){
}
L2_SH=calc_res(sectionvzdata,sectionvz,sectionvzdiff,sectionvzdiffold,ntr,ns,LNORM,L2_SH,0,1,swstestshot,ntr_glob,recpos_loc,nsrc_glob,ishot,iter,srcpos,recpos);
if(swstestshot==1){energy_SH=calc_energy(sectionvzdata,ntr,ns,energy_SH, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);}
energy_all_shots_SH=calc_energy(sectionvzdata,ntr,ns,energy_all_shots_SH, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
L2_all_shots_SH=calc_misfit(sectionvzdata,sectionvz,ntr,ns,LNORM,L2_all_shots_SH,0,1,1, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
energy_all_shots_SH=calc_energy(sectionvzdata,ntr,ns,energy_all_shots_SH, ntr_glob, recpos_loc, nsrc_glob, ishot,iter,srcpos,recpos);
}
// Tracekill
......@@ -2073,30 +2124,6 @@ int main(int argc, char **argv){
saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,sectionpdiff,sectionpdiff,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,3);
}
}
/* Write synthetic filtered seismogramms to disk */
if (SEISMO && TIME_FILT && WRITE_FILTERED_DATA==2){
if(WAVETYPE==1 || WAVETYPE==3){
if ((ADJOINT_TYPE==1)|| (ADJOINT_TYPE==3)){
catseis(sectionvx, fulldata_vx, recswitch, ntr_glob, MPI_COMM_NTR);
}
if ((ADJOINT_TYPE==1)|| (ADJOINT_TYPE==2)){
catseis(sectionvy, fulldata_vy, recswitch, ntr_glob, MPI_COMM_NTR);
}
if (ADJOINT_TYPE==4){
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_NTR);
}
}
if(WAVETYPE==2 || WAVETYPE==3){
catseis(sectionvz, fulldata_vz, recswitch, ntr_glob, MPI_COMM_NTR);
}
if(myid_ntr==0){
saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,sectionpdiff,sectionpdiff,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,4);
}
}
}
......@@ -2235,7 +2262,7 @@ int main(int argc, char **argv){
update_s_elastic_PML_SH(1, NX, 1, NY, pvz,psxz,psyz,uxz,uyz,hc,infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half,psi_vzx, psi_vzy,puipjp,pu,prho);
}
} else {
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, u, ppi, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, ppi, absorb_coeff, prho, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
}
......@@ -2299,6 +2326,7 @@ int main(int argc, char **argv){
waveconv_shot[j][i]+= (forward_prop_x[j][i][NTDTINV-hin+1]+forward_prop_y[j][i][NTDTINV-hin+1])*(psxx[j][i]+psyy[j][i]);
}else{
waveconv_shot[j][i]+= (forward_prop_p[j][i][NTDTINV-hin+1])*(psp[j][i]);
if(RTM==1) waveconv_rtm_shot[j][i]+= (forward_prop_p[j][i][NTDTINV-hin+1])*(psp[j][i]);
}
if(!ACOUSTIC){
......@@ -2384,6 +2412,7 @@ int main(int argc, char **argv){
/* interpolate unknown values */
if((IDXI>1)||(IDYI>1)){
interpol(IDXI,IDYI,waveconv_shot,1);
if(RTM==1) interpol(IDXI,IDYI,waveconv_rtm_shot,1);
}
/* calculate complete gradient */
......@@ -2392,6 +2421,8 @@ int main(int argc, char **argv){
waveconv_lam[j][i] = - DT * waveconv_shot[j][i];
if(RTM==1) waveconv_rtm_shot[j][i] = DT * waveconv_rtm_shot[j][i];
if(PARAMETERIZATION==1){
if(!ACOUSTIC)
muss = prho[j][i] * pu[j][i] * pu[j][i];
......@@ -2399,13 +2430,11 @@ int main(int argc, char **argv){
muss = 0;
lamss = prho[j][i] * ppi[j][i] * ppi[j][i] - 2.0 * muss;
if(!ACOUSTIC)
waveconv_lam[j][i] = (1.0/(4.0 * (lamss+muss) * (lamss+muss))) * waveconv_lam[j][i];
else
waveconv_lam[j][i] = (1.0/((lamss+muss) * (lamss+muss))) * waveconv_lam[j][i];
/* calculate Vp gradient */
waveconv_shot[j][i] = 2.0 * ppi[j][i] * prho[j][i] * waveconv_lam[j][i];
// waveconv_shot[j][i] = waveconv_lam[j][i];
}
if(PARAMETERIZATION==2){
......@@ -2666,6 +2695,7 @@ int main(int argc, char **argv){
if (WAVETYPE==1 || WAVETYPE==3) {
for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){
if(RTM==1) waveconv_rtm[j][i] += waveconv_rtm_shot[j][i];
waveconv[j][i] += waveconv_shot[j][i];
if(!ACOUSTIC){
waveconv_u[j][i] += waveconv_u_shot[j][i];
......@@ -2699,6 +2729,18 @@ int main(int argc, char **argv){
}
/*--------------------Reverse time migration finished-------------------------------------*/
if(RTM==1){
if (SWS_TAPER_FILE) taper_grad(waveconv_rtm,taper_coeff,srcpos,nsrc,recpos,ntr_glob,4);
sprintf(jac,"%s_rtm",JACOBIAN);
write_matrix_disk(waveconv_rtm, jac);
if(MYID==0) printf("\n Reverse time migration finished. \n\n");
MPI_Barrier(MPI_COMM_WORLD);
MPI_Buffer_detach(buff_addr,&buffsize);
MPI_Finalize();
return 0;
}
if(FORWARD_ONLY==0){
/* ----------------------------------------- */
/* Applying and output approx hessian */
......@@ -2793,7 +2835,6 @@ int main(int argc, char **argv){
if(MYID==0&&(WAVETYPE==3)) printf("\n SH: L2=%f",L2sum_all_shots_SH/energy_sum_all_shots_SH);
}
sum_killed_traces=0;
MPI_Allreduce(&killed_traces,&sum_killed_traces,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
sum_killed_traces_testshots=0;
......@@ -2804,27 +2845,6 @@ int main(int argc, char **argv){
case 2:
L2t[1]=0.0; L2t[4]=0.0;
if(JOINT_EQUAL_WEIGHTING){
if(JOINT_EQUAL_new_max){
JOINT_EQUAL_PSV=L2sum/energy_sum;
JOINT_EQUAL_SH=L2sum_SH/energy_sum_SH;
JOINT_EQUAL_PSV_all=L2sum_all_shots/energy_sum_all_shots;
JOINT_EQUAL_SH_all=L2sum_all_shots_SH/energy_sum_all_shots_SH;
JOINT_EQUAL_new_max=0;
}
L2t[1]+=(L2sum/energy_sum)/JOINT_EQUAL_PSV;
L2t[4]+=(L2sum_all_shots/energy_sum_all_shots)/JOINT_EQUAL_PSV_all;
L2t[1]+=(L2sum_SH/energy_sum_SH)/JOINT_EQUAL_SH;
L2t[4]+=(L2sum_all_shots_SH/energy_sum_all_shots_SH)/JOINT_EQUAL_SH_all;
break;
}
if(WAVETYPE==1||WAVETYPE==3){
L2t[1]+=L2sum/energy_sum;
L2t[4]+=L2sum_all_shots/energy_sum_all_shots;
......@@ -2834,7 +2854,6 @@ int main(int argc, char **argv){
L2t[1]+=L2sum_SH/energy_sum_SH;
L2t[4]+=L2sum_all_shots_SH/energy_sum_all_shots_SH;
}
if(MYID==0&&(WAVETYPE==3)) printf("\n Sum: L2=%f",L2t[4]);
break;
......@@ -3108,13 +3127,9 @@ int main(int argc, char **argv){
if(wolfe_SLS_failed) {
if (TIME_FILT==0){
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old,GAMMA);}
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old);}
else{
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f \t %f \n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old,F_LOW_PASS,GAMMA);
}
if(WAVETYPE==3 && MYID==0){
fprintf(FPL2_JOINT,"%d \t %f \t %f\n",iter,L2sum_all_shots/energy_sum_all_shots,L2sum_all_shots_SH/energy_sum_all_shots_SH);
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old,F_LOW_PASS);
}
/* No update is done here, however model fils are written to disk for easy post processing */
......@@ -3153,13 +3168,9 @@ int main(int argc, char **argv){
float diff=0.0;
diff=fabs((L2_hist[iter-2]-L2_hist[iter])/L2_hist[iter-2]);
if (TIME_FILT==0){
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f \n",alpha_SL,iter,wolfe_sum_FWI,diff,countstep-1,L2_SL_old,L2_SL_new,GAMMA);}
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \n",alpha_SL,iter,wolfe_sum_FWI,diff,countstep-1,L2_SL_old,L2_SL_new);}
else{
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f \t %f\n",alpha_SL,iter,wolfe_sum_FWI,diff,countstep-1,L2_SL_old,L2_SL_new,F_LOW_PASS,GAMMA);
}
if(WAVETYPE==3 && MYID==0){
fprintf(FPL2_JOINT,"%d \t %f \t %f\n",iter,L2sum_all_shots/energy_sum_all_shots,L2sum_all_shots_SH/energy_sum_all_shots_SH);
if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",alpha_SL,iter,wolfe_sum_FWI,diff,countstep-1,L2_SL_old,L2_SL_new,F_LOW_PASS);
}
/* initiate variables for next iteration */
......@@ -3279,25 +3290,6 @@ int main(int argc, char **argv){
fprintf(FP,"\n ***** Starting simulation (test-forward model) no. %d for shot %d of %d (rel. step length %.5f) \n",itest,ishot,nshots,eps_scale);
fprintf(FP,"\n=================================================================================================\n\n");
if (SEISMO && READREC==2){
if (ntr>0) {
dealloc_sections(ntr,ns,recpos_loc,sectionvx,sectionvy,sectionvz,sectionp,sectionpnp1,sectionpn,sectioncurl,sectiondiv,
sectionpdata,sectionpdiff,sectionpdiffold,sectionvxdata,sectionvxdiff,sectionvxdiffold,sectionvydata,
sectionvydiff,sectionvydiffold,sectionvzdata,sectionvzdiff,sectionvzdiffold);
}
free_imatrix(recpos,1,3,1,ntr_glob);
recpos=receiver(&ntr, srcpos, ishot);
recpos_loc = splitrec(recpos,&ntr_loc, ntr, recswitch);
ntr_glob=ntr;
ntr=ntr_loc;
if (ntr>0){
alloc_sections(ntr,ns,&sectionvx,&sectionvy,&sectionvz,&sectionp,&sectionpnp1,&sectionpn,&sectioncurl,&sectiondiv,
&sectionpdata,&sectionpdiff,&sectionpdiffold,&sectionvxdata,&sectionvxdiff,&sectionvxdiffold,&sectionvydata,
&sectionvydiff,&sectionvydiffold,&sectionvzdata,&sectionvzdiff,&sectionvzdiffold);
}
}
for (nt=1;nt<=8;nt++) srcpos1[nt][1]=srcpos[nt][ishot];
/*-----------------------------------*/
......@@ -3450,7 +3442,7 @@ int main(int argc, char **argv){
}
}
else /* acoustic */
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, u, ppinp1, absorb_coeff, prhonp1, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
update_p_PML(1, NX, 1, NY, pvx, pvy, psp, ppinp1, absorb_coeff, prhonp1, hc, infoout, K_x, a_x, b_x, K_x_half, a_x_half, b_x_half, K_y, a_y, b_y, K_y_half, a_y_half, b_y_half, psi_vxx, psi_vyy, psi_vxy, psi_vyx);
}
if (MYID==0){
......@@ -3620,15 +3612,6 @@ int main(int argc, char **argv){
case 2:
L2t[itest]=0.0;
if(JOINT_EQUAL_WEIGHTING){
L2t[itest]+=(L2sum/energy_sum)/JOINT_EQUAL_PSV;
L2t[itest]+=(L2sum_SH/energy_sum_SH)/JOINT_EQUAL_SH;
break;
}
if(WAVETYPE==1||WAVETYPE==3){
L2t[itest]+=L2sum/energy_sum;
}
......@@ -3825,12 +3808,9 @@ int main(int argc, char **argv){
if(MYID==0&&(!VERBOSE)){printf("L2-Norm[4] = %e (final L2 of all shots)\n",L2t[4]);}
printf("MYID = %d \t epst1[1] = %e \t epst1[2] = %e \t epst1[3] = %e \n",MYID,epst1[1],epst1[2],epst1[3]);
if (TIME_FILT==0){
fprintf(FPL2,"%e \t %e \t %e \t %e \t %e \t %e \t %e \t %e \t %f \n",opteps_vp,