diff --git a/src/IFOS2D.c b/src/IFOS2D.c index a4e55391654a55eaf87369d7e51da11bcbdc82d0..16dbe5318796ffc9efc87efe8d9677dae6945c5a 100644 --- a/src/IFOS2D.c +++ b/src/IFOS2D.c @@ -1437,6 +1437,7 @@ int main(int argc, char **argv){ if (MYID == 0) fprintf(FP,"\n start wavefield separation..."); if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1); if (MYID == 0) fprintf(FP,"\n wavefield separation finished"); + MPI_Barrier(MPI_COMM_WORLD); inseis(fprec,ishot,sectionread,ntr_glob,ns,3,iter); h=1; for(i=1;i<=ntr;i++){ @@ -1468,7 +1469,8 @@ int main(int argc, char **argv){ timedomain_filt(sectionvy_obs,FC,ORDER,ntr_glob,ns,1); } if (ADJOINT_TYPE==4){ - inseis(fprec,ishot,sectionp_obs,ntr_glob,ns,9,iter); + if (WAVESEP==1) {inseis(fprec,ishot,sectionp_obs,ntr_glob,ns,5,iter); + }else {inseis(fprec,ishot,sectionp_obs,ntr_glob,ns,9,iter);} timedomain_filt(sectionp_obs,FC,ORDER,ntr_glob,ns,1); } } @@ -1965,21 +1967,24 @@ 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); /* wavefield separation */ if (WAVESEP==1) { if (MYID == 0) fprintf(FP,"\n start wavefield separation..."); if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1); if (MYID == 0) fprintf(FP,"wavefield separation finished."); - inseis(fprec,ishot,sectionread,ntr_glob,ns,3,iter); - h=1; - for(i=1;i<=ntr;i++){ - for(j=1;j<=ns;j++){ - sectionp[h][j]=sectionread[recpos_loc[3][i]][j]; + if (FORWARD_ONLY==0) { + MPI_Barrier(MPI_COMM_WORLD); + inseis(fprec,ishot,sectionread,ntr_glob,ns,3,iter); + h=1; + for(i=1;i<=ntr;i++){ + for(j=1;j<=ns;j++){ + sectionp[h][j]=sectionread[recpos_loc[3][i]][j]; + } + h++; } - h++; - } + } } break; @@ -2056,7 +2061,8 @@ int main(int argc, char **argv){ /* read seismic data from SU file p */ /* --------------------------------- */ if(ADJOINT_TYPE==4){ /* if ADJOINT_TYPE */ - inseis(fprec,ishot,sectionread,ntr_glob,ns,9,iter); + if (WAVESEP==1) { inseis(fprec,ishot,sectionread,ntr_glob,ns,5,iter); + }else { inseis(fprec,ishot,sectionread,ntr_glob,ns,9,iter);} if ((TIME_FILT==1 )|| (TIME_FILT==2)){ timedomain_filt(sectionread,FC,ORDER,ntr_glob,ns,1); } @@ -3481,12 +3487,13 @@ int main(int argc, char **argv){ } /* wavefield separation */ if (WAVESEP==1) { - catseis(sectionvz, fulldata_vz, recswitch, ntr_glob, MPI_COMM_WORLD); + catseis(sectionvy, fulldata_vy, recswitch, ntr_glob, MPI_COMM_WORLD); catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD); if (MYID == 0) fprintf(FP,"\n start wavefield separation..."); if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,2); if (MYID == 0) fprintf(FP,"\n wavefield separation finished"); - inseis(fprec,ishot,sectionread,ntr_glob,ns,3,iter); + MPI_Barrier(MPI_COMM_WORLD); + inseis(fprec,ishot,sectionread,ntr_glob,ns,4,iter); h=1; for(i=1;i<=ntr;i++){ for(j=1;j<=ns;j++){ @@ -3555,7 +3562,8 @@ int main(int argc, char **argv){ /* read seismic data from SU file p */ /* --------------------------------- */ if(ADJOINT_TYPE==4){ /* if ADJOINT_TYPE */ - inseis(fprec,ishot,sectionread,ntr_glob,ns,9,iter); + if (WAVESEP==1) {inseis(fprec,ishot,sectionread,ntr_glob,ns,5,iter); + }else {inseis(fprec,ishot,sectionread,ntr_glob,ns,9,iter);} if ((TIME_FILT==1 )|| (TIME_FILT==2)){ timedomain_filt(sectionread,FC,ORDER,ntr_glob,ns,1); } diff --git a/src/inseis.c b/src/inseis.c index a6ddffcccf82107ef2d311c91f7456d9f4b3178a..012ae5f15b3b90a9257b6b0ee8e936054fa9b86e 100644 --- a/src/inseis.c +++ b/src/inseis.c @@ -40,15 +40,19 @@ void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int sprintf(data,"%s_vy.su.shot%d",DATA_DIR,comp); } - if(sws==3){ /* open seismic data pup */ - sprintf(data,"%s_pup.su.shot%i.it%i",SEIS_FILE,comp,iter); + if(sws==3){ /* open forward modeled seismic data pup*/ + sprintf(data,"%s_pup.su.shot%i",SEIS_FILE,comp); } - if(sws==4){ /* open seismic data pup used for step length calculation */ - sprintf(data,"%s_pup.su.step.shot%i.it%i",SEIS_FILE,comp,iter); + if(sws==4){ /* open forward modeled seismic data pup used for step length calculation */ + sprintf(data,"%s_pup.su.step.shot%i",SEIS_FILE,comp); } - /* sws 5 -- 6 not used */ + if(sws==5){ /* open seismic data pup */ + sprintf(data,"%s_pup.su.shot%i",DATA_DIR,comp); + } + + /* SWS==6 not used */ if(sws==7){ /* open convolved seismic data vx */ sprintf(data,"%s_vx.su.conv.shot%d",DATA_DIR,comp); @@ -66,7 +70,6 @@ void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int sprintf(data,"%s_vz.su.shot%d",DATA_DIR,comp); } - fpdata = fopen(data,"r"); if (fpdata==NULL) { if(MYID==0) printf(" Was not able to read %s",data); diff --git a/src/pup.c b/src/pup.c index 86bf7d0cc299e323a5f9b162e7ef6a8dd5142501..2efe945e3e465da26d773ca4f48d111eee6d31d1 100644 --- a/src/pup.c +++ b/src/pup.c @@ -51,19 +51,21 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, ny_pot = (int)(pow(2,(ceil(log(ns*2)/log(2))))); /* test: write p-data to file*/ + if (ishot==1 && iter==1) { sprintf(pf1,"%s_data_p_.shot%i.it%i.bin",SEIS_FILE,ishot, iter); file = fopen(pf1,"wb"); for (i=1;i<=ntr_glob;i++) for (j=1;j<=ns;j++) fwrite(&data_p[i][j],sizeof(float),1,file); fclose(file); - + } /* test: write vz-data to file*/ + if (ishot==1 && iter==1) { sprintf(pf1,"%s_data_vz_.shot%i.it%i.bin",SEIS_FILE,ishot, iter); file = fopen(pf1,"wb"); for (i=1;i<=ntr_glob;i++) for (j=1;j<=ns;j++) fwrite(&data_vy[i][j],sizeof(float),1,file); fclose(file); - + } /* temporary array */ filtwk=matrix(1,ny_pot/2,1,nx_pot/2); filtwk_full=matrix(1,ny_pot,1,nx_pot); @@ -145,21 +147,24 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, fft2(data_pt, data_pim, ny_pot, nx_pot,1); /* write spectrum to file */ + if (ishot==1 && iter==1) { sprintf(pf2,"%s_spectrum.p.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter); file = fopen(pf2,"wb"); for (i=1;i<=nx_pot;i++) for (j=1;j<=ny_pot;j++) fwrite(&data_pt[j][i],sizeof(float),1,file); fclose(file); + } fft2(data_vyt, data_vyim, ny_pot, nx_pot,1); /* write spectrum to file */ + if (ishot==1 && iter==1) { sprintf(pf2,"%s_spectrum.vz.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter); file = fopen(pf2,"wb"); for (i=1;i<=nx_pot;i++) for (j=1;j<=ny_pot;j++) fwrite(&data_vyt[j][i],sizeof(float),1,file); fclose(file); - + } fprintf(FP,"Perform wavefield separation... \n"); @@ -171,12 +176,13 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, } } /* write spectrum to file */ - sprintf(pf2,"%s_spectrum.pup.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter); - file = fopen(pf2,"wb"); - for (i=1;i<=nx_pot;i++) - for (j=1;j<=ny_pot;j++) fwrite(&data_pup[j][i],sizeof(float),1,file); - fclose(file); - + if (ishot==1 && iter==1) { + sprintf(pf2,"%s_spectrum.pup.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter); + file = fopen(pf2,"wb"); + for (i=1;i<=nx_pot;i++) + for (j=1;j<=ny_pot;j++) fwrite(&data_pup[j][i],sizeof(float),1,file); + fclose(file); + } fprintf(FP,"Reverse 2D FFT...\n "); @@ -186,12 +192,15 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, fprintf(FP,"Write PUP data to file...\n "); /* write spectrum to file */ - sprintf(pf2,"%s_data_pup.shot%i.it%i.bin",SEIS_FILE,ishot, iter); - file = fopen(pf2,"wb"); - for (i=1;i<=nx_pot;i++) - for (j=1;j<=ny_pot;j++) fwrite(&data_pup[j][i],sizeof(float),1,file); - fclose(file); + if (ishot==1 && iter==1) { + sprintf(pf2,"%s_data_pup.shot%i.it%i.bin",SEIS_FILE,ishot, iter); + file = fopen(pf2,"wb"); + for (i=1;i<=nx_pot;i++) + for (j=1;j<=ny_pot;j++) fwrite(&data_pup[j][i],sizeof(float),1,file); + fclose(file); + } + fprintf(FP,"Write PUP-data into p-data array...\n "); /* overwrite p data */ for (k=1;k<=ntr_glob;k++) { for (i=1;i<=ns;i++) { @@ -201,14 +210,13 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos, /* save pup*/ if (sw == 1){ - sprintf(pf,"%s_pup.su.shot%i.it%i",SEIS_FILE,ishot,iter); + sprintf(pf,"%s_pup.su.shot%i",SEIS_FILE,ishot); } if (sw == 2){ /* for step length calculation*/ - sprintf(pf,"%s_pup.su.step.shot%i.it%i",SEIS_FILE,ishot,iter); + sprintf(pf,"%s_pup.su.step.shot%i",SEIS_FILE,ishot); } outseis_glob(fp,fopen(pf,"w"), 0, data_p,recpos,recpos_loc,ntr,srcpos,1,ns,SEIS_FORMAT,ishot,1); - fprintf(FP,"Write PUP-data into p-data array...\n "); free_matrix(filtwk, 1,ny_pot/2,1,nx_pot/2); diff --git a/src/read_par_json.c b/src/read_par_json.c index 7f47dd0645048d69c941baa8354100f0f687a22d..56f7ff0349e96e1307fcb12b34d0b4d9de73c0ae 100644 --- a/src/read_par_json.c +++ b/src/read_par_json.c @@ -395,15 +395,15 @@ void read_par_json(FILE *fp, char *fileinp){ /*================================= section wavefield separation =================================*/ - if (SEISMO != 5 && SEISMO !=4 && WAVESEP == 1) { - WAVESEP=0; - fprintf(fp,"Variable SEISMO has to be set to 4 or 5 to use wavefield separation. WAVESEP is set now to %i.\n",WAVESEP); - } + if (get_int_from_objectlist("WAVESEP",number_readobjects,&WAVESEP,varname_list, value_list)){ WAVESEP=0; fprintf(fp,"Variable WAVESEP is set to default value %i.\n",WAVESEP);} else { - + if (SEISMO < 4 && WAVESEP == 1) { + WAVESEP=0; + fprintf(fp,"Variable SEISMO has to be set to 4 or 5 to use wavefield separation. WAVESEP is set now to %i.\n",WAVESEP); + } if (WAVESEP==1) { if (get_float_from_objectlist("VEL",number_readobjects,&VEL,varname_list, value_list)){ VEL=1500;