Commit a1bfa1e9 authored by laura.gassner's avatar laura.gassner
Browse files

Added option VELOCITY for acoustic inversion

parent 57a02ac2
...@@ -58,7 +58,7 @@ int main(int argc, char **argv){ ...@@ -58,7 +58,7 @@ int main(int argc, char **argv){
float L2_SH, L2sum_SH, L2_all_shots_SH, L2sum_all_shots_SH; float L2_SH, L2sum_SH, L2_all_shots_SH, L2sum_all_shots_SH;
// Pointer for dynamic wavefields: // Pointer for dynamic wavefields:
float ** psxx, ** psxy, ** psyy, ** psxz, ** psyz, **psp, ** ux, ** uy, ** uxy, ** uyx, ** Vp0, ** uttx, ** utty, ** Vs0, ** Rho0; float ** psxx, ** psxy, ** psyy, ** psxz, ** psyz, **psp, ** ux, ** uy, ** uxy, ** uyx, ** u, ** 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 ** 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_shot, **waveconv_u_shot, **waveconv_rho_shot, **waveconv_u_shot_z, **waveconv_rho_shot_z;
float ** pvxp1, ** pvyp1, ** pvzp1, ** pvxm1, ** pvym1, ** pvzm1; float ** pvxp1, ** pvyp1, ** pvzp1, ** pvxm1, ** pvym1, ** pvzm1;
...@@ -408,6 +408,8 @@ int main(int argc, char **argv){ ...@@ -408,6 +408,8 @@ int main(int argc, char **argv){
uxz = matrix(-nd+1,NY+nd,-nd+1,NX+nd); uxz = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
uyz = 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) { switch (WAVETYPE) {
...@@ -1303,7 +1305,7 @@ int main(int argc, char **argv){ ...@@ -1303,7 +1305,7 @@ int main(int argc, char **argv){
if(!ACOUSTIC) { 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); 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 { } else {
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); 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);
} }
} }
if (WAVETYPE==2 || WAVETYPE==3) { if (WAVETYPE==2 || WAVETYPE==3) {
...@@ -1746,7 +1748,7 @@ int main(int argc, char **argv){ ...@@ -1746,7 +1748,7 @@ int main(int argc, char **argv){
if(!ACOUSTIC) { 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); 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 { } else {
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); 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);
} }
} }
if (WAVETYPE==2 || WAVETYPE==3) { if (WAVETYPE==2 || WAVETYPE==3) {
...@@ -1836,7 +1838,11 @@ int main(int argc, char **argv){ ...@@ -1836,7 +1838,11 @@ int main(int argc, char **argv){
}else{ }else{
for (j=1;j<=NY;j=j+IDYI){ for (j=1;j<=NY;j=j+IDYI){
for (i=1;i<=NX;i=i+IDXI){ for (i=1;i<=NX;i=i+IDXI){
forward_prop_p[j][i][hin]=psp[j][i]; if(VELOCITY==0){
forward_prop_p[j][i][hin]=psp[j][i];
}else{
forward_prop_p[j][i][hin]=u[j][i];
}
} }
} }
} }
...@@ -2254,7 +2260,7 @@ int main(int argc, char **argv){ ...@@ -2254,7 +2260,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); 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 { } else {
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); 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);
} }
} }
...@@ -3417,7 +3423,7 @@ int main(int argc, char **argv){ ...@@ -3417,7 +3423,7 @@ int main(int argc, char **argv){
} }
} }
else /* acoustic */ else /* acoustic */
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); 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);
} }
if (MYID==0){ if (MYID==0){
......
...@@ -499,7 +499,7 @@ void update_v_acoustic_PML(int nx1, int nx2, int ny1, int ny2, int nt, ...@@ -499,7 +499,7 @@ void update_v_acoustic_PML(int nx1, int nx2, int ny1, int ny2, int nt,
float ** psi_sxx_x, float ** psi_syy_y); float ** psi_sxx_x, float ** psi_syy_y);
void update_p_PML(int nx1, int nx2, int ny1, int ny2, void update_p_PML(int nx1, int nx2, int ny1, int ny2,
float ** vx, float ** vy, float ** sp, float ** pi, float ** absorb_coeff, float **rho, float *hc, int infoout, float ** vx, float ** vy, float ** sp, float ** u, float ** pi, float ** absorb_coeff, float **rho, float *hc, int infoout,
float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half, float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half, float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,
float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx); float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx);
......
...@@ -475,9 +475,9 @@ void read_par_json(FILE *fp, char *fileinp){ ...@@ -475,9 +475,9 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_int_from_objectlist("VELOCITY",number_readobjects,&VELOCITY,varname_list, value_list)){ if (get_int_from_objectlist("VELOCITY",number_readobjects,&VELOCITY,varname_list, value_list)){
VELOCITY=0; VELOCITY=0;
fprintf(fp,"Variable VELOCITY is set to default value %d.\n",VELOCITY); fprintf(fp,"Variable VELOCITY is set to default value %d.\n",VELOCITY);
}else }/*else
if(ACOUSTIC && VELOCITY==1) if(ACOUSTIC && VELOCITY==1)
declare_error("For acoustic inversion option VELOCITY not available"); declare_error("For acoustic inversion option VELOCITY not available");*/
if (get_int_from_objectlist("USE_WORKFLOW",number_readobjects,&USE_WORKFLOW,varname_list, value_list)){ if (get_int_from_objectlist("USE_WORKFLOW",number_readobjects,&USE_WORKFLOW,varname_list, value_list)){
USE_WORKFLOW=0; USE_WORKFLOW=0;
......
...@@ -24,7 +24,7 @@ ...@@ -24,7 +24,7 @@
#include "fd.h" #include "fd.h"
void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy, float ** sp, float ** pi, float ** absorb_coeff, float **rho, float *hc, int infoout, void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy, float ** sp, float ** u, float ** pi, float ** absorb_coeff, float **rho, float *hc, int infoout,
float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half, float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half,
float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half, float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half,
float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx){ float ** psi_vxx, float ** psi_vyy, float ** psi_vxy, float ** psi_vyx){
...@@ -99,6 +99,7 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy ...@@ -99,6 +99,7 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy
sp[j][i] += g*(vxx+vyy); sp[j][i] += g*(vxx+vyy);
u[j][i] = (g/DT)*(vxx+vyy);
} }
} }
break; break;
...@@ -150,7 +151,8 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy ...@@ -150,7 +151,8 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy
g = rho[j][i] * (pi[j][i] * pi[j][i]);} g = rho[j][i] * (pi[j][i] * pi[j][i]);}
sp[j][i] += g*(vxx+vyy); sp[j][i] += g*(vxx+vyy);
u[j][i] = (g/DT)*(vxx+vyy);
} }
} }
break; break;
...@@ -205,7 +207,8 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy ...@@ -205,7 +207,8 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy
g = rho[j][i] * (pi[j][i] * pi[j][i]);} g = rho[j][i] * (pi[j][i] * pi[j][i]);}
sp[j][i] += g*(vxx+vyy); sp[j][i] += g*(vxx+vyy);
u[j][i] = (g/DT)*(vxx+vyy);
} }
} }
break; break;
...@@ -263,7 +266,8 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy ...@@ -263,7 +266,8 @@ void update_p_PML(int nx1, int nx2, int ny1, int ny2, float ** vx, float ** vy
g = rho[j][i] * (pi[j][i] * pi[j][i]);} g = rho[j][i] * (pi[j][i] * pi[j][i]);}
sp[j][i] += g*(vxx+vyy); sp[j][i] += g*(vxx+vyy);
u[j][i] = (g/DT)*(vxx+vyy);
} }
} }
break; break;
......
...@@ -119,20 +119,24 @@ void update_v_acoustic_PML(int nx1, int nx2, int ny1, int ny2, int nt, ...@@ -119,20 +119,24 @@ void update_v_acoustic_PML(int nx1, int nx2, int ny1, int ny2, int nt,
if(sw==0){ if(sw==0){
if (VELOCITY==0){ if (VELOCITY==0){
vxp1[j][i] = rip[j][i]*(sp_x)/DH; vxp1[j][i] = rip[j][i]*(sp_x)/DH;
vyp1[j][i] = rjp[j][i]*(sp_y)/DH;} vyp1[j][i] = rjp[j][i]*(sp_y)/DH;
}
else{ else{
vxp1[j][i] = rip[j][i]*(sp_x)/DH; vxp1[j][i] = rip[j][i]*(sp_x)/DH;
vyp1[j][i] = rjp[j][i]*(sp_y)/DH;} vyp1[j][i] = rjp[j][i]*(sp_y)/DH;
}
} }
if(sw==1){ if(sw==1){
if (VELOCITY==0){ if (VELOCITY==0){
vxp1[j][i] += vx[j][i]*DT; vxp1[j][i] += vx[j][i]*DT;
vyp1[j][i] += vy[j][i]*DT;} vyp1[j][i] += vy[j][i]*DT;
}
else{ else{
vxp1[j][i] = vx[j][i]; vxp1[j][i] = vx[j][i];
vyp1[j][i] = vy[j][i];} vyp1[j][i] = vy[j][i];
}
} }
vx[j][i] += DT*rip[j][i]*(sp_x)/DH; vx[j][i] += DT*rip[j][i]*(sp_x)/DH;
......
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