Skip to content

IFOS2D
IFOS2D
Commit d8de885d authored by niklas.thiel's avatar niklas.thiel
Browse files
Options

BUGFIX 

error occured in wavefield separation during steplength estimation
parent 0ca42abe
Hide whitespace changes
Inline Side-by-side
Showing with 60 additions and 41 deletions
src/IFOS2D.c
View file @ d8de885d
......@@ -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);
}
......
src/inseis.c
View file @ d8de885d
......@@ -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);
......
src/pup.c
View file @ d8de885d
......@@ -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);
......
src/read_par_json.c
View file @ d8de885d
......@@ -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;
......
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