Commit 72f5d9ef authored by Florian Wittkamp's avatar Florian Wittkamp

L-BFGS: Acoustic version

parent e85740fe
...@@ -46,12 +46,15 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float ...@@ -46,12 +46,15 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
char jac[225]; char jac[225];
FILE *FP_JAC; FILE *FP_JAC;
/*---------------------*/ /*---------------------*/
/* DEBUGGING */ /* DEBUGGING */
/*---------------------*/ /*---------------------*/
sprintf(jac,"%s_grad1_vs_it%d",JACOBIAN,iteration);
write_matrix_disk(grad_vs, jac); if(!ACOUSTIC) {
sprintf(jac,"%s_grad1_vs_it%d",JACOBIAN,iteration);
write_matrix_disk(grad_vs, jac);
}
sprintf(jac,"%s_grad1_rho_it%d",JACOBIAN,iteration); sprintf(jac,"%s_grad1_rho_it%d",JACOBIAN,iteration);
write_matrix_disk(grad_rho, jac); write_matrix_disk(grad_rho, jac);
...@@ -72,7 +75,7 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float ...@@ -72,7 +75,7 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
for (i=1;i<=NX;i++){ for (i=1;i<=NX;i++){
for (j=1;j<=NY;j++){ for (j=1;j<=NY;j++){
l++; l++;
y_LBFGS[w][l]=-grad_vs[j][i]*Vs_avg; /* VS */ if(!ACOUSTIC) y_LBFGS[w][l]=-grad_vs[j][i]*Vs_avg; /* VS */
if(NPAR_LBFGS>1) y_LBFGS[w][l+NX*NY]=-grad_rho[j][i]*rho_avg; /* RHO */ if(NPAR_LBFGS>1) y_LBFGS[w][l+NX*NY]=-grad_rho[j][i]*rho_avg; /* RHO */
if(NPAR_LBFGS>2) y_LBFGS[w][l+2*NX*NY]=-grad_vp[j][i]*Vp_avg; /* VP */ if(NPAR_LBFGS>2) y_LBFGS[w][l+2*NX*NY]=-grad_vp[j][i]*Vp_avg; /* VP */
} }
...@@ -95,8 +98,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float ...@@ -95,8 +98,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
if(w==0) w=N_LBFGS; if(w==0) w=N_LBFGS;
/* Debugging */ /* Debugging */
sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iteration,w,POS[1],POS[2]); if(!ACOUSTIC) {
FP_JAC=fopen(jac,"wb"); sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iteration,w,POS[1],POS[2]);
FP_JAC=fopen(jac,"wb");
}
if(MYID==0) printf("\n\n ------------ L-BFGS ---------------"); if(MYID==0) printf("\n\n ------------ L-BFGS ---------------");
if(MYID==0) printf("\n Start calculation L-BFGS update"); if(MYID==0) printf("\n Start calculation L-BFGS update");
...@@ -109,8 +114,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float ...@@ -109,8 +114,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
l++; l++;
/* VS */ /* VS */
y_LBFGS[w][l]+=grad_vs[j][i]*Vs_avg; /* add grad(i) to build grad(i)-grad(i-1) */ if(!ACOUSTIC){
q_LBFGS[l]=grad_vs[j][i]*Vs_avg; /* Normalisation */ y_LBFGS[w][l]+=grad_vs[j][i]*Vs_avg; /* add grad(i) to build grad(i)-grad(i-1) */
q_LBFGS[l]=grad_vs[j][i]*Vs_avg; /* Normalisation */
}
/* RHO */ /* RHO */
if(NPAR_LBFGS>1) { if(NPAR_LBFGS>1) {
...@@ -125,20 +132,22 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float ...@@ -125,20 +132,22 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
} }
/* Debugging */ /* Debugging */
fwrite(&y_LBFGS[w][l],sizeof(float),1,FP_JAC); if(!ACOUSTIC) fwrite(&y_LBFGS[w][l],sizeof(float),1,FP_JAC);
} }
} }
/*---------------------*/ /*---------------------*/
/* DEBUGGING */ /* DEBUGGING */
/*--------------------*/ /*---------------------*/
fclose(FP_JAC); if(!ACOUSTIC) {
MPI_Barrier(MPI_COMM_WORLD); fclose(FP_JAC);
sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin",JACOBIAN,iteration,w); MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(jac,3); sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin",JACOBIAN,iteration,w);
MPI_Barrier(MPI_COMM_WORLD); if (MYID==0) mergemod(jac,3);
sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iteration,w,POS[1],POS[2]); MPI_Barrier(MPI_COMM_WORLD);
remove(jac); sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iteration,w,POS[1],POS[2]);
remove(jac);
}
/*----------------------------------*/ /*----------------------------------*/
/* call L-BFGS Algorithm */ /* call L-BFGS Algorithm */
...@@ -158,8 +167,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float ...@@ -158,8 +167,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
l++; l++;
/* VS */ /* VS */
y_LBFGS[w][l]=-grad_vs[j][i]*Vs_avg; /* add -grad(i-1) to build grad(i)-grad(i-1) */ if(!ACOUSTIC) {
grad_vs[j][i]=r_LBFGS[l]*Vs_avg; /* Denormalization */ y_LBFGS[w][l]=-grad_vs[j][i]*Vs_avg; /* add -grad(i-1) to build grad(i)-grad(i-1) */
grad_vs[j][i]=r_LBFGS[l]*Vs_avg; /* Denormalization */
}
/* RHO */ /* RHO */
if(NPAR_LBFGS>1) { if(NPAR_LBFGS>1) {
...@@ -184,9 +195,11 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float ...@@ -184,9 +195,11 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
/*---------------------*/ /*---------------------*/
/* DEBUGGING */ /* DEBUGGING */
/*---------------------*/ /*---------------------*/
sprintf(jac,"%s_grad2_vs_it%d",JACOBIAN,iteration); if(!ACOUSTIC){
write_matrix_disk(grad_vs, jac); sprintf(jac,"%s_grad2_vs_it%d",JACOBIAN,iteration);
write_matrix_disk(grad_vs, jac);
}
sprintf(jac,"%s_grad2_rho_it%d",JACOBIAN,iteration); sprintf(jac,"%s_grad2_rho_it%d",JACOBIAN,iteration);
write_matrix_disk(grad_rho, jac); write_matrix_disk(grad_rho, jac);
...@@ -306,7 +319,8 @@ void lbfgs_reset(int iter, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS1, float ...@@ -306,7 +319,8 @@ void lbfgs_reset(int iter, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS1, float
extern int NX,NY,MYID; extern int NX,NY,MYID;
if(MYID==0) printf("\n\n ------------ L-BFGS ---------------"); if(MYID==0) printf("\n\n ------------ L-BFGS ---------------");
if(MYID==0) printf("\n Reset L-BFGS at Iteration %d",iter); if(MYID==0&&iter>1) printf("\n Reset L-BFGS at iteration %d",iter);
if(MYID==0&&iter==1) printf("\n L-BFGS will be used from iteration %d on",iter);
for(l=1;l<=N_LBFGS;l++){ for(l=1;l<=N_LBFGS;l++){
for(m=1;m<=NPAR_LBFGS*NX*NY;m++){ for(m=1;m<=NPAR_LBFGS*NX*NY;m++){
......
...@@ -98,7 +98,7 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST ...@@ -98,7 +98,7 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST
*FC=workflow[WORKFLOW_STAGE][7]; *FC=workflow[WORKFLOW_STAGE][7];
} }
} else { } else {
if(MYID==0)printf("\n TIME_FILT cannot be activated due to it is not activated in the JSON File \n"); if(MYID==0&&(workflow[WORKFLOW_STAGE][6]>0))printf("\n TIME_FILT cannot be activated due to it is not activated in the JSON File \n");
} }
/* Change of wavetype */ /* Change of wavetype */
if(wavetype_start!=3&&(WAVETYPE!=workflow[WORKFLOW_STAGE][8])){ if(wavetype_start!=3&&(WAVETYPE!=workflow[WORKFLOW_STAGE][8])){
......
...@@ -61,8 +61,10 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float * ...@@ -61,8 +61,10 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float *
w=iter%N_LBFGS; w=iter%N_LBFGS;
if(w==0) w=N_LBFGS; if(w==0) w=N_LBFGS;
sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]); if(!ACOUSTIC){
FP_JAC=fopen(jac,"wb"); sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]);
FP_JAC=fopen(jac,"wb");
}
if(LBFGS_NPAR>1){ if(LBFGS_NPAR>1){
sprintf(jac2,"%s_s_LBFGS_rho_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]); sprintf(jac2,"%s_s_LBFGS_rho_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]);
...@@ -313,9 +315,11 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float * ...@@ -313,9 +315,11 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float *
if(itest==0){ if(itest==0){
if(GRAD_METHOD==2){ if(GRAD_METHOD==2){
l++; l++;
if(!ACOUSTIC) s_LBFGS[w][l]=(unp1[j][i]-u[j][i])/Vs_avg; if(!ACOUSTIC) {
fwrite(&s_LBFGS[w][l],sizeof(float),1,FP_JAC); s_LBFGS[w][l]=(unp1[j][i]-u[j][i])/Vs_avg;
fwrite(&s_LBFGS[w][l],sizeof(float),1,FP_JAC);
}
if(LBFGS_NPAR>1) { if(LBFGS_NPAR>1) {
s_LBFGS[w][l+NX*NY]=(rhonp1[j][i]-rho[j][i])/rho_avg; s_LBFGS[w][l+NX*NY]=(rhonp1[j][i]-rho[j][i])/rho_avg;
...@@ -326,6 +330,7 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float * ...@@ -326,6 +330,7 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float *
s_LBFGS[w][l+2*NX*NY]=(pinp1[j][i]-pi[j][i])/Vp_avg; s_LBFGS[w][l+2*NX*NY]=(pinp1[j][i]-pi[j][i])/Vp_avg;
fwrite(&s_LBFGS[w][l+2*NY*NX],sizeof(float),1,FP_JAC3); fwrite(&s_LBFGS[w][l+2*NY*NX],sizeof(float),1,FP_JAC3);
} }
} }
if(!ACOUSTIC) u[j][i] = unp1[j][i]; if(!ACOUSTIC) u[j][i] = unp1[j][i];
rho[j][i] = rhonp1[j][i]; rho[j][i] = rhonp1[j][i];
...@@ -336,13 +341,15 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float * ...@@ -336,13 +341,15 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float *
} }
if(GRAD_METHOD==2&&(itest==0)){ if(GRAD_METHOD==2&&(itest==0)){
fclose(FP_JAC); if(!ACOUSTIC){
MPI_Barrier(MPI_COMM_WORLD); fclose(FP_JAC);
sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin",JACOBIAN,iter+1,w); MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(jac,3); sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin",JACOBIAN,iter+1,w);
MPI_Barrier(MPI_COMM_WORLD); if (MYID==0) mergemod(jac,3);
sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]); MPI_Barrier(MPI_COMM_WORLD);
remove(jac); sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]);
remove(jac);
}
if(LBFGS_NPAR>1){ if(LBFGS_NPAR>1){
fclose(FP_JAC2); fclose(FP_JAC2);
......
...@@ -2887,25 +2887,31 @@ int main(int argc, char **argv){ ...@@ -2887,25 +2887,31 @@ int main(int argc, char **argv){
if(GRAD_FILTER==1){smooth(waveconv,1,1,Vs_avg,FC);} if(GRAD_FILTER==1){smooth(waveconv,1,1,Vs_avg,FC);}
} }
if (SWS_TAPER_FILE){ /* read taper from BIN-File*/ if (SWS_TAPER_FILE && !ACOUSTIC){ /* read taper from BIN-File*/
taper_grad(waveconv_u,taper_coeff,srcpos,nsrc,recpos,ntr_glob,5);} taper_grad(waveconv_u,taper_coeff,srcpos,nsrc,recpos,ntr_glob,5);}
if (SWS_TAPER_FILE){ /* read taper from BIN-File*/ if (SWS_TAPER_FILE){ /* read taper from BIN-File*/
taper_grad(waveconv_rho,taper_coeff,srcpos,nsrc,recpos,ntr_glob,6);} taper_grad(waveconv_rho,taper_coeff,srcpos,nsrc,recpos,ntr_glob,6);}
if(GRAD_FILTER==1){smooth(waveconv_u,2,1,Vs_avg,FC);} if(GRAD_FILTER==1&& !ACOUSTIC){smooth(waveconv_u,2,1,Vs_avg,FC);}
if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,FC);} if(GRAD_FILTER==1){smooth(waveconv_rho,3,1,Vs_avg,FC);}
if(WOLFE_CONDITION) { if(WOLFE_CONDITION) {
for (j=1;j<=NY;j=j+IDY){ for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){ for (i=1;i<=NX;i=i+IDX){
if(WAVETYPE!=2) waveconv_old[j][i]=waveconv[j][i];
waveconv_u_old[j][i] = waveconv_u[j][i];
waveconv_rho_old[j][i] = waveconv_rho[j][i];
if(WAVETYPE!=2){
waveconv_old[j][i]=waveconv[j][i];
ppinp1[j][i] = ppi[j][i];
}
if(!ACOUSTIC){
waveconv_u_old[j][i] = waveconv_u[j][i];
punp1[j][i] =pu[j][i];
}
waveconv_rho_old[j][i] = waveconv_rho[j][i];
prhonp1[j][i] = prho[j][i]; prhonp1[j][i] = prho[j][i];
if(WAVETYPE!=2) ppinp1[j][i] = ppi[j][i];
punp1[j][i] =pu[j][i];
} }
} }
} }
...@@ -2921,7 +2927,7 @@ int main(int argc, char **argv){ ...@@ -2921,7 +2927,7 @@ int main(int argc, char **argv){
for (j=1;j<=NY;j=j+IDY){ for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){ for (i=1;i<=NX;i=i+IDX){
if(WAVETYPE!=2) waveconv_up[j][i]=waveconv[j][i]; if(WAVETYPE!=2) waveconv_up[j][i]=waveconv[j][i];
waveconv_u_up[j][i] = waveconv_u[j][i]; if(!ACOUSTIC) waveconv_u_up[j][i] = waveconv_u[j][i];
waveconv_rho_up[j][i] = waveconv_rho[j][i]; waveconv_rho_up[j][i] = waveconv_rho[j][i];
} }
} }
...@@ -3018,7 +3024,7 @@ int main(int argc, char **argv){ ...@@ -3018,7 +3024,7 @@ int main(int argc, char **argv){
for (i=1;i<=NX;i=i+IDX){ for (i=1;i<=NX;i=i+IDX){
prho[j][i]=prhonp1[j][i]; prho[j][i]=prhonp1[j][i];
if(WAVETYPE!=2) ppi[j][i]=ppinp1[j][i]; if(WAVETYPE!=2) ppi[j][i]=ppinp1[j][i];
pu[j][i]=punp1[j][i]; if(!ACOUSTIC) pu[j][i]=punp1[j][i];
} }
} }
......
...@@ -701,16 +701,11 @@ void read_par_json(FILE *fp, char *fileinp){ ...@@ -701,16 +701,11 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_int_from_objectlist("GRAD_METHOD",number_readobjects,&GRAD_METHOD,varname_list, value_list)) if (get_int_from_objectlist("GRAD_METHOD",number_readobjects,&GRAD_METHOD,varname_list, value_list))
err("Variable GRAD_METHOD could not be retrieved from the json input file!"); err("Variable GRAD_METHOD could not be retrieved from the json input file!");
else { else {
if(ACOUSTIC){
GRAD_METHOD=1;
fprintf(fp,"For acoustic modelling currently only GRAD_METHOD=%d possible.\n",GRAD_METHOD);}
if(GRAD_METHOD==1) { if(GRAD_METHOD==1) {
WOLFE_CONDITION=0; WOLFE_CONDITION=0;
} }
if(GRAD_METHOD==2) { if(GRAD_METHOD==2) {
if(ACOUSTIC){err("Currently no L-BFGS with acoustic FWI.");}
if (get_int_from_objectlist("LBFGS_SURFACE",number_readobjects,&LBFGS_SURFACE,varname_list, value_list)){ if (get_int_from_objectlist("LBFGS_SURFACE",number_readobjects,&LBFGS_SURFACE,varname_list, value_list)){
LBFGS_SURFACE=0; LBFGS_SURFACE=0;
fprintf(fp,"Variable LBFGS_SURFACE is set to default value %d.\n",LBFGS_SURFACE); fprintf(fp,"Variable LBFGS_SURFACE is set to default value %d.\n",LBFGS_SURFACE);
......
...@@ -27,6 +27,7 @@ int check_wolfe(float steplength, float misfit_old, float misfit_new, float ** g ...@@ -27,6 +27,7 @@ int check_wolfe(float steplength, float misfit_old, float misfit_new, float ** g
/* global */ /* global */
extern int NX,NY,MYID; extern int NX,NY,MYID;
extern FILE *FP; extern FILE *FP;
extern int ACOUSTIC;
/* local */ /* local */
float left=0.0, right=0.0; float left=0.0, right=0.0;
...@@ -34,11 +35,13 @@ int check_wolfe(float steplength, float misfit_old, float misfit_new, float ** g ...@@ -34,11 +35,13 @@ int check_wolfe(float steplength, float misfit_old, float misfit_new, float ** g
float temp; float temp;
float temp2; float temp2;
temp=matrix_product(grad_old_vs, update_vs); if(!ACOUSTIC){
grad_dot_p_old+=temp; temp=matrix_product(grad_old_vs, update_vs);
temp2=global_maximum(grad_old_vs); grad_dot_p_old+=temp;
if(fabs(temp2)>0) grad_dot_p_old_norm+=temp/temp2; temp2=global_maximum(grad_old_vs);
grad_dot_p_new+=matrix_product(grad_new_vs, update_vs); if(fabs(temp2)>0) grad_dot_p_old_norm+=temp/temp2;
grad_dot_p_new+=matrix_product(grad_new_vs, update_vs);
}
if(NPAR_LBFGS>1) { if(NPAR_LBFGS>1) {
temp=matrix_product(grad_old_rho, update_rho); temp=matrix_product(grad_old_rho, update_rho);
...@@ -87,7 +90,7 @@ int check_wolfe(float steplength, float misfit_old, float misfit_new, float ** g ...@@ -87,7 +90,7 @@ int check_wolfe(float steplength, float misfit_old, float misfit_new, float ** g
fprintf(FP,"\n Condition 2 NOT OK: (%.1f>= %.1f)", left, right); fprintf(FP,"\n Condition 2 NOT OK: (%.1f>= %.1f)", left, right);
return 2; return 2;
} }
fprintf(FP,"\n Steplength %.3f satisfy wolfe condition",steplength); fprintf(FP,"\n Steplength %.3f satisfy wolfe condition",steplength);
return 0; return 0;
} }
......
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