...
 
Commits (4)
......@@ -489,8 +489,11 @@ Default values are:
With the use of a workflow file, you can define different FWI stages. For instance, one FWI stage could refer to one corner frequency of a low-pass filter or to a specific time window. Every line in the workflow file corresponds to one FWI stage. To use a workflow file the switch USE\_WORKFLOW have to be set to 1. The algorithm will automatically change to the next line of the workflow file if the abort criterium of the current line is reached or if no step length could be found which reduces the misfit. The structure of the variables inside the workflow file is as follow: \\
\noindent
\begin{tabular}{llllllllllIlll}
\# & INV\_VS & INV\_VP & INV\_RHO & PRO & TIME\_FILT & F\_HIGH\_PASS & F\_LOW\_PASS & WAVETYPE & 0 & 0 & EPRECOND & EPSILON\_WE & GAMMA \\
\begin{tabular}{lllllllll}
\# & INV\_VS & INV\_VP & INV\_RHO & PRO & TIME\_FILT & F\_HIGH\_PASS & F\_LOW\_PASS & WAVETYPE\\
\end{tabular}
\begin{tabular}{llllll}
\# & 0 & 0 & EPRECOND & EPSILON\_WE & GAMMA \\
\end{tabular}
\ \\
......
......@@ -967,6 +967,7 @@ void read_par_json(FILE *fp, char *fileinp){
if (FORWARD_ONLY>0){ /* only forward modeling is applied */
FORWARD_ONLY=1;
ITERMAX=1;
INV_STF=0;
strcpy(INV_MODELFILE,MFILE);
DTINV=1;
INVTYPE=2;
......
......@@ -143,17 +143,8 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
}
if(FORWARD_ONLY==1){
if(ADJOINT_TYPE==1) {
sprintf(vxf,"%s_vx.su.shot%d",SEIS_FILE,ishot);
sprintf(vyf,"%s_vy.su.shot%d",SEIS_FILE,ishot);
}
if(ADJOINT_TYPE==2){
sprintf(vyf,"%s_vy.su.shot%d",SEIS_FILE,ishot);
}
if(ADJOINT_TYPE==3){
sprintf(vxf,"%s_vx.su.shot%d",SEIS_FILE,ishot);
}
if(WAVETYPE==2 || WAVETYPE==3) {
sprintf(vzf,"%s_vz.su.shot%d",SEIS_FILE,ishot);
}
......@@ -165,7 +156,7 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
switch (SEISMO){
case 1 : /* particle velocities only */
if (WAVETYPE==1 || WAVETYPE==3) {
if(ADJOINT_TYPE==1) {
if(ADJOINT_TYPE==1 || ADJOINT_TYPE==0) {
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
......@@ -201,7 +192,7 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
case 4 : /* everything */
if (WAVETYPE==1 || WAVETYPE==3) {
if(ADJOINT_TYPE==1) {
if(ADJOINT_TYPE==1 || ADJOINT_TYPE==0) {
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf);
......@@ -232,7 +223,7 @@ void saveseis_glob(FILE *fp, float **sectionvx, float **sectionvy,float **sectio
case 5 : /* everything except curl and div */
if (WAVETYPE==1 || WAVETYPE==3) {
if(ADJOINT_TYPE==1) {
if(ADJOINT_TYPE==1 || ADJOINT_TYPE==0) {
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vx) to\n\t %s \n",0,ntr,vxf);
outseis_glob(fp,fopen(vxf,"w"),1,sectionvx,recpos,recpos_loc,ntr,srcpos,nsrc,ns,SEIS_FORMAT,ishot,1);
if(VERBOSE==1) fprintf(fp," PE %d is writing %d seismograms (vy) to\n\t %s \n",0,ntr,vyf);
......
......@@ -34,7 +34,7 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx){
int i,j, m, fdoh, h, h1, l;
float vxx, vyy, vxy, vyx;
float vxx, vyy;
float dhi, dthalbe;
extern float DT, DH;
extern int MYID, FDORDER, FW, L;
......@@ -62,21 +62,85 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
for (j=ny1;j<=ny2;j++){
for (i=nx1;i<=nx2;i++){
vxx = ( hc[1]*(vx[j][i] -vx[j][i-1]))*dhi;
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i]))*dhi;
vyx = ( hc[1]*(vy[j][i+1]-vy[j][i]))*dhi;
/* left boundary */
if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){
vxy = ( hc[1]*(vx[j+1][i]-vx[j][i]))*dhi;
psi_vxx[j][i] = b_x[i] * psi_vxx[j][i] + a_x[i] * vxx;
vxx = vxx / K_x[i] + psi_vxx[j][i];
}
vyy = ( hc[1]*(vy[j][i] -vy[j-1][i]))*dhi;
/* 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];
}
/* computing sums of the old memory variables */
sump=0.0;
for (l=1;l<=L;l++){
sump+=p[j][i][l];
}
/* updating components of the stress tensor, partially */
sp[j][i] += (g[j][i]*(vxx+vyy))+(dthalbe*sump);
// u[j][i] = ((g[j][i]/DT)*(vxx+vyy))-((2.0*f[j][i]/DT)*vyy)+(0.5*sump);
/* now updating the memory-variables and sum them up*/
sump=0.0;
for (l=1;l<=L;l++){
p[j][i][l] = bjm[l]*(p[j][i][l]*cjm[l]-(e[j][i][l]*(vxx+vyy)));
sump += p[j][i][l];
}
/* and now the components of the stress tensor are
completely updated */
sp[j][i]+=(dthalbe*sump);
// u[j][i]+=(0.5*sump);
}
}
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];
psi_vyx[j][i] = b_x_half[i] * psi_vyx[j][i] + a_x_half[i] * vyx;
vyx = vyx / K_x_half[i] + psi_vyx[j][i];
}
/* right boundary */
......@@ -87,22 +151,90 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
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];
}
/* computing sums of the old memory variables */
sump=0.0;
for (l=1;l<=L;l++){
sump+=p[j][i][l];
}
/* updating components of the stress tensor, partially */
sp[j][i] += (g[j][i]*(vxx+vyy))+(dthalbe*sump);
// u[j][i] = ((g[j][i]/DT)*(vxx+vyy))-((2.0*f[j][i]/DT)*vyy)+(0.5*sump);
/* now updating the memory-variables and sum them up*/
sump=0.0;
for (l=1;l<=L;l++){
p[j][i][l] = bjm[l]*(p[j][i][l]*cjm[l]-(e[j][i][l]*(vxx+vyy)));
sump += p[j][i][l];
}
/* and now the components of the stress tensor are
completely updated */
sp[j][i]+=(dthalbe*sump);
// u[j][i]+=(0.5*sump);
}
}
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)){
/*psi_vyx[j][h] = b_x_half[h] * psi_vyx[j][h] + a_x_half[h] * vyx;
vyx = vyx / K_x_half[h] + psi_vyx[j][h];*/
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];
psi_vyx[j][h1] = b_x_half[h1] * psi_vyx[j][h1] + a_x_half[h1] * vyx;
vyx = vyx / K_x_half[h1] + psi_vyx[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;
psi_vxy[j][i] = b_y_half[j] * psi_vxy[j][i] + a_y_half[j] * vxy;
vyy = vyy / K_y[j] + psi_vyy[j][i];
vxy = vxy / K_y_half[j] + psi_vxy[j][i];
}
/* bottom boundary */
......@@ -113,12 +245,82 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
psi_vyy[h1][i] = b_y[h1] * psi_vyy[h1][i] + a_y[h1] * vyy;
vyy = vyy / K_y[h1] + psi_vyy[h1][i];
}
/* computing sums of the old memory variables */
sump=0.0;
for (l=1;l<=L;l++){
sump+=p[j][i][l];
}
/* updating components of the stress tensor, partially */
sp[j][i] += (g[j][i]*(vxx+vyy))+(dthalbe*sump);
// u[j][i] = ((g[j][i]/DT)*(vxx+vyy))-((2.0*f[j][i]/DT)*vyy)+(0.5*sump);
/* now updating the memory-variables and sum them up*/
sump=0.0;
for (l=1;l<=L;l++){
p[j][i][l] = bjm[l]*(p[j][i][l]*cjm[l]-(e[j][i][l]*(vxx+vyy)));
sump += p[j][i][l];
}
/* and now the components of the stress tensor are
completely updated */
sp[j][i]+=(dthalbe*sump);
// u[j][i]+=(0.5*sump);
}
}
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_vxy[j][i] = b_y_half[j] * psi_vxy[j][i] + a_y_half[j] * vxy;
vxy = vxy / K_y_half[j] + psi_vxy[j][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];
}
psi_vxy[h1][i] = b_y_half[h1] * psi_vxy[h1][i] + a_y_half[h1] * vxy;
vxy = vxy / K_y_half[h1] + psi_vxy[h1][i];
/* 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];
}
/* computing sums of the old memory variables */
......@@ -130,6 +332,9 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
/* updating components of the stress tensor, partially */
sp[j][i] += (g[j][i]*(vxx+vyy))+(dthalbe*sump);
// u[j][i] = ((g[j][i]/DT)*(vxx+vyy))-((2.0*f[j][i]/DT)*vyy)+(0.5*sump);
/* now updating the memory-variables and sum them up*/
sump=0.0;
for (l=1;l<=L;l++){
......@@ -137,10 +342,12 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
sump += p[j][i][l];
}
/* and now the components of the stress tensor are
completely updated */
sp[j][i]+=(dthalbe*sump);
// u[j][i]+=(0.5*sump);
}
}
break;
......@@ -150,18 +357,12 @@ void update_p_visc_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float **
for (i=nx1;i<=nx2;i++){
vxx = 0.0;
vyy = 0.0;
vyx = 0.0;
vxy = 0.0;
for (m=1; m<=fdoh; m++) {
vxx += hc[m]*(vx[j][i+m-1] -vx[j][i-m] );
vyy += hc[m]*(vy[j+m-1][i] -vy[j-m][i] );
vyx += hc[m]*(vy[j][i+m] -vy[j][i-m+1]);
vxy += hc[m]*(vx[j+m][i] -vx[j-m+1][i]);
}
vxx *= dhi;
vyy *= dhi;
vyx *= dhi;
vxy *= dhi;
sump=0.0;
for (l=1;l<=L;l++){
......