Commit 6761a267 authored by laura.gassner's avatar laura.gassner

minor addaptions, removal of erroneous option for STF inversion

parent 3363b15f
......@@ -32,6 +32,7 @@ void model_acoustic(float **rho, float **pi){
extern char MFILE[STRING_SIZE];
extern char INV_MODELFILE[STRING_SIZE];
extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
extern int SWS_TAPER_FILE;
/* local variables */
float piv, vp, rhov, taperv, **taper;
......@@ -54,81 +55,80 @@ void model_acoustic(float **rho, float **pi){
flvp=vector(1,nodes);
fltaper=vector(1,nodes);
taper=matrix(-nd+1,NY+nd,-nd+1,NX+nd);
if(SWS_TAPER_FILE) taper=matrix(-nd+1,NY+nd,-nd+1,NX+nd);
flfile=fopen("model_true/flnodes.toy_example_ac.start","r");
if (flfile==NULL) err(" FL-file could not be opened !");
/* Read parameters */
for (l=1;l<=nodes;l++){
fgets(cline,255,flfile);
if (cline[0]!='#'){
sscanf(cline,"%f%f%f%f",&fldepth[l], &flrho[l], &flvp[l], &fltaper[l]);
}
else l=l-1;
}
/* Print parameters */
if(MYID==0){
printf(" ------------------------------------------------------------------ \n\n");
printf(" Information of FL nodes: \n\n");
printf(" \t depth \t vp \n\n");
for (l=1;l<=nodes;l++){
printf(" \t %f \t %f\n\n",fldepth[l],flvp[l]);
}
printf(" ------------------------------------------------------------------ \n\n");
printf(" ------------------------------------------------------------------ \n\n");
printf(" Information of FL nodes: \n\n");
printf(" \t depth \t vp \n\n");
for (l=1;l<=nodes;l++){
printf(" \t %f \t %f\n\n",fldepth[l],flvp[l]);
}
printf(" ------------------------------------------------------------------ \n\n");
}
/* loop over global grid */
for (i=1;i<=NXG;i++){
for (l=1;l<nodes;l++){
if (fldepth[l]==fldepth[l+1]){
if ((i==1) && (MYID==0)){
printf("depth: %f m: double node\n",fldepth[l]);}}
else{
for (j=(int)(fldepth[l]/DH)+1;j<=(int)(fldepth[l+1]/DH);j++){
vp=0.0;
vp=(DH*(j-1)-fldepth[l])*(flvp[l+1]-flvp[l])/(fldepth[l+1]-fldepth[l])+flvp[l];
vp=vp*1000.0;
piv=vp;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
/* loop over global grid - vp */
for (i=1;i<=NXG;i++){
for (l=1;l<nodes;l++){
if (fldepth[l]==fldepth[l+1]){
if ((i==1) && (MYID==0)){
printf("depth: %f m: double node\n",fldepth[l]);
}
}else{
for (j=(int)(fldepth[l]/DH)+1;j<=(int)(fldepth[l+1]/DH);j++){
vp=0.0;
vp=(DH*(j-1)-fldepth[l])*(flvp[l+1]-flvp[l])/(fldepth[l+1]-fldepth[l])+flvp[l];
vp=vp*1000.0;
piv=vp;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
pi[jj][ii]=piv;
}
}
}
}
}
for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
vp=0.0;
vp=flvp[nodes]*1000.0;
piv=vp;
vp=0.0;
vp=flvp[nodes]*1000.0;
piv=vp;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
pi[jj][ii]=piv;
}
pi[jj][ii]=piv;
}
}
}
/* loop over global grid - taper */
if(SWS_TAPER_FILE){
for (i=1;i<=NXG;i++){
for (l=1;l<nodes;l++){
if (fldepth[l]==fldepth[l+1]){
......@@ -155,7 +155,7 @@ void model_acoustic(float **rho, float **pi){
}
for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
taperv=0.0;
taperv=fltaper[nodes];
......@@ -170,6 +170,7 @@ void model_acoustic(float **rho, float **pi){
}
}
}
}
free_vector(fldepth,1,nodes);
free_vector(flrho,1,nodes);
......@@ -178,7 +179,7 @@ void model_acoustic(float **rho, float **pi){
/**************************************************/
/* creation of density models */
/* creation of density models from true model */
/**************************************************/
/*read FL nodes from File*/
......@@ -189,6 +190,7 @@ void model_acoustic(float **rho, float **pi){
flfile=fopen("model_true/flnodes.toy_example_ac","r");
if (flfile==NULL) err(" FL-file could not be opened !");
/* Read parameters */
for (l=1;l<=nodes;l++){
fgets(cline,255,flfile);
if (cline[0]!='#'){
......@@ -197,19 +199,19 @@ void model_acoustic(float **rho, float **pi){
else l=l-1;
}
/* Print parameters */
if(MYID==0){
printf(" ------------------------------------------------------------------ \n\n");
printf(" Information of FL nodes: \n\n");
printf(" \t depth \t vp \t rho \n\n");
for (l=1;l<=nodes;l++){
printf(" \t %f \t %f \t %f \n\n",fldepth[l],flvp[l],flrho[l]);
}
printf(" ------------------------------------------------------------------ \n\n");
printf(" ------------------------------------------------------------------ \n\n");
printf(" Information of FL nodes: \n\n");
printf(" \t depth \t vp \t rho \n\n");
for (l=1;l<=nodes;l++){
printf(" \t %f \t %f \t %f \n\n",fldepth[l],flvp[l],flrho[l]);
}
printf(" ------------------------------------------------------------------ \n\n");
}
/* loop over global grid */
/* loop over global grid - density */
for (i=1;i<=NXG;i++){
for (l=1;l<nodes;l++){
if (fldepth[l]==fldepth[l+1]){
......@@ -275,22 +277,24 @@ void model_acoustic(float **rho, float **pi){
sprintf(modfile,"%s_vp_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s",TAPER_FILE_NAME);
writemod(modfile,taper,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s.%i%i",TAPER_FILE_NAME,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s",TAPER_FILE_NAME_RHO);
writemod(modfile,taper,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s.%i%i",TAPER_FILE_NAME_RHO,POS[1],POS[2]);
remove(modfile);
if(SWS_TAPER_FILE){
sprintf(modfile,"%s",TAPER_FILE_NAME);
writemod(modfile,taper,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s.%i%i",TAPER_FILE_NAME,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s",TAPER_FILE_NAME_RHO);
writemod(modfile,taper,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s.%i%i",TAPER_FILE_NAME_RHO,POS[1],POS[2]);
remove(modfile);
}
free_matrix(taper,-nd+1,NY+nd,-nd+1,NX+nd);
if(SWS_TAPER_FILE) free_matrix(taper,-nd+1,NY+nd,-nd+1,NX+nd);
}
......@@ -36,7 +36,7 @@
"RUN_MULTIPLE_SHOTS" : "1",
"Acoustic Computation" : "comment",
"ACOUSTIC" : "1",
"ACOUSTIC" : "0",
"Model" : "comment",
"READMOD" : "0",
......
......@@ -36,7 +36,7 @@
"RUN_MULTIPLE_SHOTS" : "1",
"Acoustic Computation" : "comment",
"ACOUSTIC" : "1",
"ACOUSTIC" : "0",
"Model" : "comment",
"READMOD" : "0",
......
......@@ -43,4 +43,4 @@ make clean
rm jacobian/toy_example/*toy_example*.*.*.*
rm model/*.bin.*.*
rm model/toy_example/*.bin.*.*
rm taper/*.bin.*.*
\ No newline at end of file
rm *.bin.*.*
\ No newline at end of file
......@@ -71,7 +71,7 @@ int main(int argc, char **argv){
** sectionvy_conv=NULL, ** sectionvy_obs=NULL, ** sectionvx_conv=NULL,** sectionvx_obs=NULL, ** sectionp_conv=NULL,** sectionp_obs=NULL, * source_time_function=NULL;
float ** absorb_coeff, ** taper_coeff, * epst1, * epst2, * epst3, * picked_times;
float ** srcpos=NULL, **srcpos_loc=NULL, ** srcpos1=NULL, **srcpos_loc_back=NULL, ** signals=NULL, ** signals_rec=NULL, *hc=NULL, ** dsignals=NULL;
float ** srcpos=NULL, **srcpos_loc=NULL, ** srcpos1=NULL, **srcpos_loc_back=NULL, ** signals=NULL, ** signals_rec=NULL, *hc=NULL;
int ** recpos=NULL, ** recpos_loc=NULL;
int * DTINV_help;
......@@ -700,9 +700,6 @@ int main(int argc, char **argv){
case 2: FC_EXT=filter_frequencies(&nfrq); FC=FC_EXT[FREQ_NR]; break;
}
if (INV_STF==1)
dsignals=fmatrix(1,nsrc,1,NT);
QUELLART_OLD = QUELLART;
for(iter=1;iter<=ITERMAX;iter++){ /* fullwaveform iteration loop */
......@@ -946,23 +943,7 @@ int main(int argc, char **argv){
signals=NULL;
signals=wavelet(srcpos_loc,nsrc_loc,ishot);
if((iter>1)&&(nsrc_loc>0)&&(INV_STF==1)){
/* find maximum of dsignals */
sig_max = 0.0;
for(iq=1;iq<=NT;iq++){
if(fabs(dsignals[ishot][iq])>sig_max){
sig_max = fabs(dsignals[ishot][iq]);
}
}
/* update wavelet */
for(iq=1;iq<=NT;iq++){
signals[1][iq] -= 0.05 * (dsignals[ishot][iq]/sig_max);
}
}
/* initialize wavefield with zero */
if (L){
if(!ACOUSTIC)
......@@ -994,7 +975,6 @@ int main(int argc, char **argv){
for (nt=1;nt<=NT;nt++){
/* Check if simulation is still stable */
/*if (isnan(pvy[NY/2][NX/2])) err(" Simulation is unstable !");*/
if (isnan(pvy[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !");}
......@@ -1353,25 +1333,8 @@ int main(int argc, char **argv){
}
MPI_Barrier(MPI_COMM_WORLD);
if((iter>1)&&(nsrc_loc>0)&&(INV_STF==1)){
/* find maximum of dsignals */
sig_max = 0.0;
for(iq=1;iq<=NT;iq++){
if(fabs(dsignals[ishot][iq])>sig_max){
sig_max = fabs(dsignals[ishot][iq]);
}
}
/* update wavelet */
for(iq=1;iq<=NT;iq++){
signals[1][iq] -= 0.05 * (dsignals[ishot][iq]/sig_max);
}
}
/* initialize wavefield with zero */
if (L){
if(!ACOUSTIC)
......@@ -1444,7 +1407,6 @@ int main(int argc, char **argv){
for (nt=1;nt<=NT;nt++){
/* Check if simulation is still stable */
/*if (isnan(pvy[NY/2][NX/2])) err(" Simulation is unstable !");*/
if (isnan(pvy[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !");
......@@ -1872,7 +1834,6 @@ int main(int argc, char **argv){
for (nt=1;nt<=NT;nt++){
/* Check if simulation is still stable */
/*if (isnan(pvy[NY/2][NX/2])) err(" Simulation is unstable !");*/
if (isnan(pvy[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !");
......@@ -1963,18 +1924,6 @@ int main(int argc, char **argv){
if (infoout) fprintf(FP," finished (real time: %4.2f s).\n",time7-time6);
}
/* calculate change of the source wavelet */
if ((nsrc_loc>0)&&(INV_STF==1)){
for (lq=1;lq<=nsrc_loc;lq++) {
iq=(int)srcpos_loc[1][lq];
jq=(int)srcpos_loc[2][lq];
if(!ACOUSTIC)
dsignals[ishot][invtimer] = (-psxx[jq][lq]-psyy[jq][lq])/2.0;
else
dsignals[ishot][invtimer] = (-psp[jq][lq])/2.0;
}
}
/*if(nt==hin1){*/
if(DTINV_help[NT-nt+1]==1){
......@@ -2587,7 +2536,9 @@ int main(int argc, char **argv){
for (nt=1;nt<=NT;nt++){
/* Check if simulation is still stable */
if (isnan(pvy[NY/2][NX/2])) err(" Simulation is unstable !");
if (isnan(pvy[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !");}
infoout = !(nt%10000);
......@@ -3351,7 +3302,6 @@ int main(int argc, char **argv){
if (nsrc_loc>0){
free_matrix(signals,1,nsrc_loc,1,NT);
if(INV_STF==1) free_matrix(dsignals,1,nsrc_loc,1,NT);
free_matrix(srcpos_loc,1,8,1,nsrc_loc);
free_matrix(srcpos_loc_back,1,6,1,nsrc_loc);
}
......
......@@ -29,26 +29,17 @@ void inseis_source_wavelet(float *section, int ns, int ishot){
extern int MYID;
extern char SIGNAL_FILE[STRING_SIZE];
/* declaration of local variables */
int j;
float dump;
segy tr;
char data[STRING_SIZE];
FILE *fpdata;
extern float TIME, DH, DT, REFREC[4];
sprintf(data,"%s.shot%d",SIGNAL_FILE,ishot);
fpdata = fopen(data,"r");
/* declaration of local variables */
int i,j;
segy tr;
float dump;
if (fpdata==NULL) err(" Source wavelet not found ");
/* SEGY (without file-header) */
fread(&tr,240,1,fpdata);
......
......@@ -243,7 +243,7 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
conv_FD(wavelet,source_time_function,stf_conv_wavelet,ns);
if(TAPER_STF)
taper(stf_conv_wavelet, ns, FC);
taper(stf_conv_wavelet, ns, fc);
// /* --------------- writing out the observed seismograms --------------- */
// sprintf(obs_y_tmp,"%s.shot%d.obs",SEIS_FILE_P,ishot);
......
......@@ -34,8 +34,8 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy
int i,j, m, fdoh, h, h1;
float g;
float vxx, vyy, vxy, vyx;
float dhi;
float vxx, vyy;
float dhi;
extern float DT, DH;
extern int MYID, FDORDER, INVMAT1, FW;
extern int FREE_SURF, BOUNDARY;
......@@ -105,6 +105,171 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy
}
}
break;
case 4:
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vxx = ( hc[1]*(vx[j][i] -vx[j][i-1])
+ hc[2]*(vx[j][i+1]-vx[j][i-2]))*dhi;
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i])
+ hc[2]*(vy[j+1][i]-vy[j-2][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
psi_vxx[j][i] = b_x[i] * psi_vxx[j][i] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i];
}
/* right boundary */
if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){
h1 = (i-nx2+2*FW);
h = i;
psi_vxx[j][h1] = b_x[h1] * psi_vxx[j][h1] + a_x[h1] * vxx;
vxx = vxx / K_x[h1] + psi_vxx[j][h1];
}
/* top boundary */
if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){
psi_vyy[j][i] = b_y[j] * psi_vyy[j][i] + a_y[j] * vyy;
vyy = vyy / K_y[j] + psi_vyy[j][i];
}
/* bottom boundary */
if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){
h1 = (j-ny2+2*FW);
h = j;
psi_vyy[h1][i] = b_y[h1] * psi_vyy[h1][i] + a_y[h1] * vyy;
vyy = vyy / K_y[h1] + psi_vyy[h1][i];
}
if (INVMAT1==1){
g = rho[j][i] * (pi[j][i] * pi[j][i]);}
sp[j][i] += g*(vxx+vyy);
}
}
break;
case 6:
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vxx = ( hc[1]*(vx[j][i] -vx[j][i-1])
+ hc[2]*(vx[j][i+1]-vx[j][i-2])
+ hc[3]*(vx[j][i+2]-vx[j][i-3]))*dhi;
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i])
+ hc[2]*(vy[j+1][i]-vy[j-2][i])
+ hc[3]*(vy[j+2][i]-vy[j-3][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
psi_vxx[j][i] = b_x[i] * psi_vxx[j][i] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i];
}
/* right boundary */
if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){
h1 = (i-nx2+2*FW);
h = i;
psi_vxx[j][h1] = b_x[h1] * psi_vxx[j][h1] + a_x[h1] * vxx;
vxx = vxx / K_x[h1] + psi_vxx[j][h1];
}
/* top boundary */
if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){
psi_vyy[j][i] = b_y[j] * psi_vyy[j][i] + a_y[j] * vyy;
vyy = vyy / K_y[j] + psi_vyy[j][i];
}
/* bottom boundary */
if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){
h1 = (j-ny2+2*FW);
h = j;
psi_vyy[h1][i] = b_y[h1] * psi_vyy[h1][i] + a_y[h1] * vyy;
vyy = vyy / K_y[h1] + psi_vyy[h1][i];
}
if (INVMAT1==1){
g = rho[j][i] * (pi[j][i] * pi[j][i]);}
sp[j][i] += g*(vxx+vyy);
}
}
break;
case 8:
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vxx = ( hc[1]*(vx[j][i] -vx[j][i-1])
+ hc[2]*(vx[j][i+1]-vx[j][i-2])
+ hc[3]*(vx[j][i+2]-vx[j][i-3])
+ hc[4]*(vx[j][i+3]-vx[j][i-4]))*dhi;
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i])
+ hc[2]*(vy[j+1][i]-vy[j-2][i])
+ hc[3]*(vy[j+2][i]-vy[j-3][i])
+ hc[4]*(vy[j+3][i]-vy[j-4][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
psi_vxx[j][i] = b_x[i] * psi_vxx[j][i] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i];
}
/* right boundary */
if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){
h1 = (i-nx2+2*FW);
h = i;
psi_vxx[j][h1] = b_x[h1] * psi_vxx[j][h1] + a_x[h1] * vxx;
vxx = vxx / K_x[h1] + psi_vxx[j][h1];
}
/* top boundary */
if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){
psi_vyy[j][i] = b_y[j] * psi_vyy[j][i] + a_y[j] * vyy;
vyy = vyy / K_y[j] + psi_vyy[j][i];
}
/* bottom boundary */
if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){
h1 = (j-ny2+2*FW);
h = j;
psi_vyy[h1][i] = b_y[h1] * psi_vyy[h1][i] + a_y[h1] * vyy;
vyy = vyy / K_y[h1] + psi_vyy[h1][i];
}
if (INVMAT1==1){
g = rho[j][i] * (pi[j][i] * pi[j][i]);}
sp[j][i] += g*(vxx+vyy);
}
}
break;
default:
for (j=ny1;j<=ny2;j++){
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
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