From 72f5d9ef11c49f2ec34630d008ab00082afaa711 Mon Sep 17 00:00:00 2001 From: Florian Wittkamp Date: Fri, 4 Dec 2015 13:36:23 +0100 Subject: [PATCH] L-BFGS: Acoustic version --- src/LBFGS.c | 60 +++++++++++++++++++++++--------------- src/apply_workflow.c | 2 +- src/calc_mat_change_test.c | 31 ++++++++++++-------- src/denise.c | 24 +++++++++------ src/read_par_json.c | 7 +---- src/wolfe_condition.c | 15 ++++++---- 6 files changed, 82 insertions(+), 57 deletions(-) diff --git a/src/LBFGS.c b/src/LBFGS.c index b433ef8..0431411 100644 --- a/src/LBFGS.c +++ b/src/LBFGS.c @@ -46,12 +46,15 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float char jac[225]; FILE *FP_JAC; - + /*---------------------*/ /* 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); 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 for (i=1;i<=NX;i++){ for (j=1;j<=NY;j++){ 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>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 if(w==0) w=N_LBFGS; /* Debugging */ - 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(!ACOUSTIC) { + 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 Start calculation L-BFGS update"); @@ -109,8 +114,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float l++; /* VS */ - 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 */ + if(!ACOUSTIC){ + 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 */ if(NPAR_LBFGS>1) { @@ -125,20 +132,22 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float } /* Debugging */ - fwrite(&y_LBFGS[w][l],sizeof(float),1,FP_JAC); + if(!ACOUSTIC) fwrite(&y_LBFGS[w][l],sizeof(float),1,FP_JAC); } } /*---------------------*/ /* DEBUGGING */ - /*--------------------*/ - fclose(FP_JAC); - MPI_Barrier(MPI_COMM_WORLD); - sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin",JACOBIAN,iteration,w); - if (MYID==0) mergemod(jac,3); - MPI_Barrier(MPI_COMM_WORLD); - sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iteration,w,POS[1],POS[2]); - remove(jac); + /*---------------------*/ + if(!ACOUSTIC) { + fclose(FP_JAC); + MPI_Barrier(MPI_COMM_WORLD); + sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin",JACOBIAN,iteration,w); + if (MYID==0) mergemod(jac,3); + MPI_Barrier(MPI_COMM_WORLD); + 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 */ @@ -158,8 +167,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float l++; /* VS */ - 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 */ + if(!ACOUSTIC) { + 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 */ if(NPAR_LBFGS>1) { @@ -184,9 +195,11 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float /*---------------------*/ /* DEBUGGING */ /*---------------------*/ - sprintf(jac,"%s_grad2_vs_it%d",JACOBIAN,iteration); - write_matrix_disk(grad_vs, jac); - + if(!ACOUSTIC){ + sprintf(jac,"%s_grad2_vs_it%d",JACOBIAN,iteration); + write_matrix_disk(grad_vs, jac); + } + sprintf(jac,"%s_grad2_rho_it%d",JACOBIAN,iteration); 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 extern int NX,NY,MYID; 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(m=1;m<=NPAR_LBFGS*NX*NY;m++){ diff --git a/src/apply_workflow.c b/src/apply_workflow.c index 2939049..88e309b 100644 --- a/src/apply_workflow.c +++ b/src/apply_workflow.c @@ -98,7 +98,7 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST *FC=workflow[WORKFLOW_STAGE][7]; } } 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 */ if(wavetype_start!=3&&(WAVETYPE!=workflow[WORKFLOW_STAGE][8])){ diff --git a/src/calc_mat_change_test.c b/src/calc_mat_change_test.c index 13c1cf7..1222035 100644 --- a/src/calc_mat_change_test.c +++ b/src/calc_mat_change_test.c @@ -61,8 +61,10 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float * w=iter%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]); - FP_JAC=fopen(jac,"wb"); + if(!ACOUSTIC){ + 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){ 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 * if(itest==0){ if(GRAD_METHOD==2){ l++; - - if(!ACOUSTIC) s_LBFGS[w][l]=(unp1[j][i]-u[j][i])/Vs_avg; - fwrite(&s_LBFGS[w][l],sizeof(float),1,FP_JAC); + + if(!ACOUSTIC) { + 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) { 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 * 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); } + } if(!ACOUSTIC) u[j][i] = unp1[j][i]; rho[j][i] = rhonp1[j][i]; @@ -336,13 +341,15 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float * } if(GRAD_METHOD==2&&(itest==0)){ - fclose(FP_JAC); - MPI_Barrier(MPI_COMM_WORLD); - sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin",JACOBIAN,iter+1,w); - if (MYID==0) mergemod(jac,3); - MPI_Barrier(MPI_COMM_WORLD); - sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]); - remove(jac); + if(!ACOUSTIC){ + fclose(FP_JAC); + MPI_Barrier(MPI_COMM_WORLD); + sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin",JACOBIAN,iter+1,w); + if (MYID==0) mergemod(jac,3); + MPI_Barrier(MPI_COMM_WORLD); + 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){ fclose(FP_JAC2); diff --git a/src/denise.c b/src/denise.c index 26e3751..7944ee2 100644 --- a/src/denise.c +++ b/src/denise.c @@ -2887,25 +2887,31 @@ int main(int argc, char **argv){ 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);} if (SWS_TAPER_FILE){ /* read taper from BIN-File*/ 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(WOLFE_CONDITION) { for (j=1;j<=NY;j=j+IDY){ 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]; - 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){ for (j=1;j<=NY;j=j+IDY){ for (i=1;i<=NX;i=i+IDX){ 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]; } } @@ -3018,7 +3024,7 @@ int main(int argc, char **argv){ for (i=1;i<=NX;i=i+IDX){ prho[j][i]=prhonp1[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]; } } diff --git a/src/read_par_json.c b/src/read_par_json.c index 3d05e50..92aec0a 100644 --- a/src/read_par_json.c +++ b/src/read_par_json.c @@ -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)) err("Variable GRAD_METHOD could not be retrieved from the json input file!"); else { - if(ACOUSTIC){ - GRAD_METHOD=1; - fprintf(fp,"For acoustic modelling currently only GRAD_METHOD=%d possible.\n",GRAD_METHOD);} if(GRAD_METHOD==1) { WOLFE_CONDITION=0; } 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)){ LBFGS_SURFACE=0; fprintf(fp,"Variable LBFGS_SURFACE is set to default value %d.\n",LBFGS_SURFACE); diff --git a/src/wolfe_condition.c b/src/wolfe_condition.c index 1821781..d33f35c 100644 --- a/src/wolfe_condition.c +++ b/src/wolfe_condition.c @@ -27,6 +27,7 @@ int check_wolfe(float steplength, float misfit_old, float misfit_new, float ** g /* global */ extern int NX,NY,MYID; extern FILE *FP; + extern int ACOUSTIC; /* local */ float left=0.0, right=0.0; @@ -34,11 +35,13 @@ int check_wolfe(float steplength, float misfit_old, float misfit_new, float ** g float temp; float temp2; - temp=matrix_product(grad_old_vs, update_vs); - grad_dot_p_old+=temp; - temp2=global_maximum(grad_old_vs); - if(fabs(temp2)>0) grad_dot_p_old_norm+=temp/temp2; - grad_dot_p_new+=matrix_product(grad_new_vs, update_vs); + if(!ACOUSTIC){ + temp=matrix_product(grad_old_vs, update_vs); + grad_dot_p_old+=temp; + temp2=global_maximum(grad_old_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) { 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 fprintf(FP,"\n Condition 2 NOT OK: (%.1f>= %.1f)", left, right); return 2; } - + fprintf(FP,"\n Steplength %.3f satisfy wolfe condition",steplength); return 0; } -- 2.26.2