diff --git a/doc/latex/6_Parameter_Definition.tex b/doc/latex/6_Parameter_Definition.tex index 1880ac56d16d438ad99009aa3941186ce5219414..2f425eb09532ddb9316c39d92d8c1102eb138463 100644 --- a/doc/latex/6_Parameter_Definition.tex +++ b/doc/latex/6_Parameter_Definition.tex @@ -535,6 +535,8 @@ The corresponding functions are copied out of DENISE Black Edition, which is mai "WOLFE_CONDITION" : "1", "WOLFE_TRY_OLD_STEPLENGTH" : "1", "WOLFE_NUM_TEST" : "5", + "WOLFE_C1_SL" : "1e-4", + "WOLFE_C2_SL" : "0.9", "LBFGS_STEP_LENGTH" : "0", \end{verbatim}}} @@ -561,7 +563,7 @@ The ensure a stable convergence when L-BFGS is used the step length has to satis \bigtriangledown f(x+\alpha \cdot p)^{T} \cdot p \geq c_2 \cdot\bigtriangledown f(x)^{T} \cdot p \end{align} To verify this condition the full gradient (with all shots!) of the updated model have to be calculated, so the verification of this condition is very expensive in terms of computation time. In theory always a step length of one should be tried first, however to save time it is possible to chose the step length from the iteration before which satisfied the wolfe condition. This behavior can be controlled with WOLFE\_TRY\_OLD\_STEPLENGTH, if this is set to 1 always the old step length will be tried first. In practice the experience has shown that at the beginning of a FWI stage a couple of step lengths will be tested and then the algorithm converges to a step length which will be accepted by the wolfe condition until the misfit cannot be reduced further. With the variable WOLFE\_NUM\_TEST a maximum number of test calculations can be defined. If after WOLFE\_NUM\_TEST tests no step length could be found, the algorithm takes the step length of the tests which reduced the misfit most, even if this step length does not satisfy the wolfe condition. If no step length was found which reduces the misfit, the algorithm will switch to the next FWI stage in the workflow file or abort the inversion. - The variable $c_1$ will be $10^{-7}\cdot max(f(x)^{-1}$ and $c_2$ is set to 0.9 (they are hard coded). + The variable $c_1$ will be $10^{-7}\cdot max(f(x)^{-1}$ and $c_2$ is set to 0.9. With WOLFE\_C1\_SL and WOLFE\_C2\_SL it is possible to manipulate $c_1$ and $c_2$ with the JSON input file. If the parameters are not set, they will be set to default values, so in most cases it should be the best to do not specify theses parameters in the JSON file. A step length search which is based on the wolfe condition can be used which the switch WOLFE\_CONDITION. Please note that this step length search is only available if GRAD\_METHOD==2. diff --git a/par/in_and_out/DENISE_INV_all_parameters.json b/par/in_and_out/DENISE_INV_all_parameters.json index 5b3f3af3ab34345b27bca30765f57dfaa9e201fa..8f47e2d9e9ced4fc6dc0ef3bd9aa9c3f276eefde 100644 --- a/par/in_and_out/DENISE_INV_all_parameters.json +++ b/par/in_and_out/DENISE_INV_all_parameters.json @@ -111,12 +111,15 @@ "WOLFE_CONDITION" : "1", "WOLFE_NUM_TEST" : "5", "WOLFE_TRY_OLD_STEPLENGTH" : "0", + "WOLFE_C1_SL" : "1e-4", + "WOLFE_C2_SL" : "0.9", "Approx. Hessian" : "comment", "EPRECOND" : "3", "EPSILON_WE" : "0.005", "EPRECOND_ITER" : "0", + "EPRECOND_PER_SHOT" : "1", "Workflow" : "comment", "USE_WORKFLOW" : "0", diff --git a/src/LBFGS.c b/src/LBFGS.c index a0b121bac0a80deb93d2e39b618a80f56fbf3607..df97134650d975cabfaa8e472996e592c8632faa 100644 --- a/src/LBFGS.c +++ b/src/LBFGS.c @@ -245,7 +245,7 @@ void lbfgs_core(int iteration, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS, flo /* give output so stdout */ if(VERBOSE || 1) { for(w=1;w<=N_LBFGS;w++) { - fprintf(FP,"\n rho_LBFGS(%d)=%e\n",w,rho_LBFGS[w]); + fprintf(FP,"\n rho_LBFGS(%2d)=%e",w,rho_LBFGS[w]); } fprintf(FP,"\n h0=%e\n",h0); } diff --git a/src/denise.c b/src/denise.c index de0a8f4aea3926f631b88b5d289ae39966b4d759..b889f63264352803b0932f86042fc40aa87459c4 100644 --- a/src/denise.c +++ b/src/denise.c @@ -111,7 +111,7 @@ int main(int argc, char **argv){ int LBFGS=0,LBFGS_NPAR=3; int LBFGS_iter_start=1; float LBFGS_L2_temp; - float **s_LBFGS,**bfgsmod2,**bfgsmod3,**y_LBFGS, **bfgsgrad2, **bfgsgrad3, *rho_LBFGS,*bfgsscale2,*bfgsscale3; + float **s_LBFGS,**y_LBFGS, *rho_LBFGS; int l=0; int w=0; int m=0; @@ -581,6 +581,10 @@ int main(int argc, char **argv){ waveconv_rho_shot = matrix(-nd+1,NY+nd,-nd+1,NX+nd); if(WOLFE_CONDITION){ + + c1_SL=WOLFE_C1_SL; + c2_SL=WOLFE_C2_SL; + waveconv_old= matrix(-nd+1,NY+nd,-nd+1,NX+nd); if(!ACOUSTIC) waveconv_u_old= matrix(-nd+1,NY+nd,-nd+1,NX+nd); waveconv_rho_old= matrix(-nd+1,NY+nd,-nd+1,NX+nd); @@ -1095,6 +1099,14 @@ int main(int argc, char **argv){ if(MYID==0){ if(iter==1){ FPL2=fopen(MISFIT_LOG_FILE,"w"); + /* Write header for misfit log file */ + if(GRAD_METHOD==1) { + if (TIME_FILT==0){ + fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \n");} + else{ + fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t FC \n"); + } + } } if(iter>1){FPL2=fopen(MISFIT_LOG_FILE,"a");} } @@ -1518,7 +1530,7 @@ int main(int argc, char **argv){ fprintf(FP,"\n==================================================================================\n"); fprintf(FP,"\n MYID=%d * Starting simulation (forward model) for shot %d of %d. Iteration %d ** \n",MYID,ishot,nshots,iter); - fprintf(FP,"\n==================================================================================\n\n"); + fprintf(FP,"\n==================================================================================\n"); for (nt=1;nt<=8;nt++) srcpos1[nt][1]=srcpos[nt][ishot]; @@ -1561,13 +1573,12 @@ int main(int argc, char **argv){ /*------------------------------------------------------------------------------*/ if (((TIME_FILT==1) || (TIME_FILT==2)) && (QUELLART!=6) && (INV_STF==0)){ - fprintf(FP," Time Domain Filter applied: lowpass filter, corner frequency of %.2f Hz, order %d\n",FC,ORDER); + fprintf(FP,"\n Time Domain Filter applied: Lowpass with corner frequency of %.2f Hz, order %d\n",FC,ORDER); /*time domain filtering of the source signal */ if(WAVETYPE==1||WAVETYPE==3) timedomain_filt(signals,FC,ORDER,nsrc_loc,ns,1); if(WAVETYPE==2||WAVETYPE==3) timedomain_filt(signals_SH,FC,ORDER,nsrc_loc,ns,1); - if(WAVETYPE==1||WAVETYPE==3){ if ((QUELLTYPB==1)|| (QUELLTYPB==2)){ /*time domain filtering of the observed data sectionvy_obs */ @@ -2122,7 +2133,7 @@ int main(int argc, char **argv){ if(MYID==0){ printf("\n==================================================================================\n"); printf("\n MYID=%d ***** Starting simulation (backward model) for shot %d of %d ********** \n",MYID,irec,nshots1); - printf("\n==================================================================================\n\n"); + printf("\n==================================================================================\n"); } @@ -2972,6 +2983,8 @@ int main(int argc, char **argv){ use_wolfe_failsafe=1; break; } else { + fprintf(FP,"\n After %d simulations no step length could be found which reduces the misfit and satisfy the wolfe condition.",countstep); + fprintf(FP,"\n Will continue without model update."); wolfe_SLS_failed=1; break; } @@ -3028,6 +3041,21 @@ int main(int argc, char **argv){ if(INVMAT!=0) break; } + if(wolfe_SLS_failed) { + + if (TIME_FILT==0){ + if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old);} + else{ + if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old,FC); + } + + /* No update is done here, however model fils are written to disk for easy post processing */ + alpha_SL=0.0; + calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,INVMAT,alpha_SL,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start); + + alpha_SL_old=1; + } + /*-----------------------------------------------------*/ /* Wolfe condition: Model update */ /*-----------------------------------------------------*/ @@ -3082,6 +3110,7 @@ int main(int argc, char **argv){ gradient_optimization=1; } + opteps_vp=0.0; opteps_vs=0.0; opteps_rho=0.0; @@ -3828,11 +3857,12 @@ int main(int argc, char **argv){ break; } - if(MYID==0) printf("Changing to corner frequency of %4.2f Hz \n",FC); + if(MYID==0) printf("\n Changing to corner frequency of %4.2f Hz \n",FC); /* Restart L-BFGS at next iteration */ LBFGS_iter_start=iter+1; wolfe_SLS_failed=0; + alpha_SL_old=1; } /* ------------------------------------------------- */ @@ -3858,11 +3888,12 @@ int main(int argc, char **argv){ s=1; min_iter_help=0; min_iter_help=iter+MIN_ITER; - if(MYID==0) printf("Changing to corner frequency of %4.2f Hz \n",FC); + if(MYID==0) printf("\n Changing to corner frequency of %4.2f Hz \n",FC); /* Restart L-BFGS at next iteration */ LBFGS_iter_start=iter+1; wolfe_SLS_failed=0; + alpha_SL_old=1; } } diff --git a/src/exchange_par.c b/src/exchange_par.c index 75e9b60305a60f3bb29464864b44514666d55982..8289159bd6614502788b8e33f4c65325945dc3c0 100644 --- a/src/exchange_par.c +++ b/src/exchange_par.c @@ -110,6 +110,8 @@ void exchange_par(void){ extern int WOLFE_CONDITION; extern int WOLFE_NUM_TEST; extern int WOLFE_TRY_OLD_STEPLENGTH; + extern float WOLFE_C1_SL; + extern float WOLFE_C2_SL; /* definition of local variables */ @@ -200,6 +202,9 @@ void exchange_par(void){ fdum[62]=EPSILON_WE; + fdum[63]=WOLFE_C1_SL; + fdum[64]=WOLFE_C2_SL; + idum[1] = NPROCX; idum[2] = NPROCY; @@ -465,6 +470,9 @@ void exchange_par(void){ EPSILON_WE=fdum[62]; + WOLFE_C1_SL=fdum[63]; + WOLFE_C2_SL=fdum[64]; + NPROCX = idum[1]; NPROCY = idum[2]; NPROC = idum[4]; diff --git a/src/globvar.h b/src/globvar.h index 93920bdc12f19a8c2b631816cbd7e39c302896ce..b0deb08be394d9e3296668702f858e8e8a62e8a5 100644 --- a/src/globvar.h +++ b/src/globvar.h @@ -96,6 +96,8 @@ int LBFGS_STEP_LENGTH; int WOLFE_CONDITION=0; int WOLFE_NUM_TEST=5; int WOLFE_TRY_OLD_STEPLENGTH=1; +float WOLFE_C1_SL; +float WOLFE_C2_SL; /* trace kill variables */ int TRKILL; diff --git a/src/read_par_json.c b/src/read_par_json.c index e6b23f2ece38a515a1d8b95164f65a80ba6365c2..7d55b0b22ef8427bad89088910dcb4ae0d0778fc 100644 --- a/src/read_par_json.c +++ b/src/read_par_json.c @@ -115,6 +115,9 @@ void read_par_json(FILE *fp, char *fileinp){ extern int WOLFE_CONDITION; extern int WOLFE_NUM_TEST; extern int WOLFE_TRY_OLD_STEPLENGTH; + extern float WOLFE_C1_SL; + extern float WOLFE_C2_SL; + /* definition of local variables */ int number_readobjects=0,fserr=0; @@ -731,6 +734,14 @@ void read_par_json(FILE *fp, char *fileinp){ WOLFE_TRY_OLD_STEPLENGTH=1; fprintf(fp,"Variable WOLFE_TRY_OLD_STEPLENGTH is set to default value %d.\n",WOLFE_TRY_OLD_STEPLENGTH); } + if (get_float_from_objectlist("WOLFE_C1_SL",number_readobjects,&WOLFE_C1_SL,varname_list, value_list)){ + WOLFE_C1_SL=1e-4; + fprintf(fp,"Variable WOLFE_C1_SL is set to default value %f.\n",WOLFE_C1_SL); + } + if (get_float_from_objectlist("WOLFE_C2_SL",number_readobjects,&WOLFE_C2_SL,varname_list, value_list)){ + WOLFE_C2_SL=0.9; + fprintf(fp,"Variable WOLFE_C2_SL is set to default value %f.\n",WOLFE_C2_SL); + } } } } diff --git a/src/wolfe_condition.c b/src/wolfe_condition.c index d33f35c54de2166f3e84c037eafae1aaf536671f..a88dd8237568a7e549a11dba26e620fd9000c536 100644 --- a/src/wolfe_condition.c +++ b/src/wolfe_condition.c @@ -59,7 +59,7 @@ int check_wolfe(float steplength, float misfit_old, float misfit_new, float ** g grad_dot_p_new+=matrix_product(grad_new_vp, update_vp); } - fprintf(FP,"\n\n--------- Wolfe condition ----------- "); + fprintf(FP,"\n\n ------- Wolfe condition -------- "); /*---------------------------------*/ /* Check condition 1 */