Skip to content

GitLab

Commit e85740fe authored by Florian Wittkamp's avatar Florian Wittkamp
Browse files

L-BFGS: Added debugging output & renamed variables 

Added debugging output for VP and renamed some variables to be consistent within the code.
parent 7797c6b2
Hide whitespace changes
Inline Side-by-side
Showing with 73 additions and 46 deletions
src/LBFGS.c
View file @ e85740fe
......@@ -34,6 +34,7 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
extern int POS[3];
extern char JACOBIAN[STRING_SIZE];
extern int WAVETYPE;
extern int ACOUSTIC;
/* local */
int m=0,v=0,w=0;
......@@ -49,12 +50,17 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
/*---------------------*/
/* DEBUGGING */
/*---------------------*/
sprintf(jac,"%s_grad1_it%d",JACOBIAN,iteration);
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);
if(WAVETYPE==1||WAVETYPE==3) {
sprintf(jac,"%s_grad1_vp_it%d",JACOBIAN,iteration);
write_matrix_disk(grad_vp, jac);
}
/*-------------------------------------------------*/
/* Init L-BFGS at iter==LBFGS_iter_start */
/*-------------------------------------------------*/
......@@ -89,7 +95,7 @@ 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_it%d_w%d.bin.%i.%i",JACOBIAN,iteration,w,POS[1],POS[2]);
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 ---------------");
......@@ -128,10 +134,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
/*--------------------*/
fclose(FP_JAC);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_y_LBFGS_it%d_w%d.bin",JACOBIAN,iteration,w);
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_it%d_w%d.bin.%i.%i",JACOBIAN,iteration,w,POS[1],POS[2]);
sprintf(jac,"%s_y_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iteration,w,POS[1],POS[2]);
remove(jac);
/*----------------------------------*/
......@@ -178,12 +184,16 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
/*---------------------*/
/* DEBUGGING */
/*---------------------*/
sprintf(jac,"%s_grad2_it%d",JACOBIAN,iteration);
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);
if(WAVETYPE==1||WAVETYPE==3) {
sprintf(jac,"%s_grad2_vp_it%d",JACOBIAN,iteration);
write_matrix_disk(grad_vp, jac);
}
}
void lbfgs_core(int iteration, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS, float ** y_LBFGS, float * rho_LBFGS,float *q_LBFGS,float *alpha_LBFGS,float *r_LBFGS) {
......
src/calc_mat_change_test.c
View file @ e85740fe
......@@ -26,7 +26,7 @@
#include "fd.h"
void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float ** waveconv_u, float ** rho, float ** rhonp1, float ** pi, float ** pinp1, float ** u, float ** unp1, int iter,
int epstest, int INVMAT, float eps_scale, int itest, int nfstart, float ** u_start, float ** pi_start, float ** rho_start,int wavetype_start,float **bfgsmod,int bfgsnum,int bfgspar,float Vs_avg,float Vp_avg,float rho_avg,int LBFGS_iter_start){
int epstest, int INVMAT, float eps_scale, int itest, int nfstart, float ** u_start, float ** pi_start, float ** rho_start,int wavetype_start,float **s_LBFGS,int N_LBFGS,int LBFGS_NPAR,float Vs_avg,float Vp_avg,float rho_avg,int LBFGS_iter_start){
/*--------------------------------------------------------------------------*/
......@@ -55,20 +55,25 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float *
extern char JACOBIAN[STRING_SIZE];
char jac[225],jac2[225];
FILE *FP_JAC,*FP_JAC2;
FILE *FP_JAC,*FP_JAC2,*FP_JAC3;
if(GRAD_METHOD==2&&(itest==0)){
w=iter%bfgsnum;
if(w==0) w=bfgsnum;
w=iter%N_LBFGS;
if(w==0) w=N_LBFGS;
sprintf(jac,"%s_bfgsmod_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]);
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(bfgspar>1){
sprintf(jac2,"%s_bfgsmod_rho_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]);
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]);
FP_JAC2=fopen(jac2,"wb");
}
if(LBFGS_NPAR>2){
sprintf(jac2,"%s_s_LBFGS_vp_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]);
FP_JAC3=fopen(jac2,"wb");
}
if(MYID==0) printf("\n\n ------------ L-BFGS ---------------");
if(MYID==0) printf("\n Saving model difference for L-BFGS");
if(MYID==0) printf("\n At Iteration %i in L-BFGS vector %i\n",iter,w);
......@@ -309,16 +314,17 @@ void calc_mat_change_test(float ** waveconv, float ** waveconv_rho, float *
if(GRAD_METHOD==2){
l++;
if(!ACOUSTIC) bfgsmod[w][l]=(unp1[j][i]-u[j][i])/Vs_avg;
fwrite(&bfgsmod[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(bfgspar>1) {
bfgsmod[w][l+NX*NY]=(rhonp1[j][i]-rho[j][i])/rho_avg;
fwrite(&bfgsmod[w][l+NY*NX],sizeof(float),1,FP_JAC2);
if(LBFGS_NPAR>1) {
s_LBFGS[w][l+NX*NY]=(rhonp1[j][i]-rho[j][i])/rho_avg;
fwrite(&s_LBFGS[w][l+NY*NX],sizeof(float),1,FP_JAC2);
}
if(bfgspar>2){
bfgsmod[w][l+2*NX*NY]=(pinp1[j][i]-pi[j][i])/Vp_avg;
if(LBFGS_NPAR>2){
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];
......@@ -332,19 +338,29 @@ 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_bfgsmod_it%d_w%d.bin",JACOBIAN,iter+1,w);
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_bfgsmod_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]);
sprintf(jac,"%s_s_LBFGS_vs_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]);
remove(jac);
if(bfgspar>1){
if(LBFGS_NPAR>1){
fclose(FP_JAC2);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_bfgsmod_rho_it%d_w%d.bin",JACOBIAN,iter+1,w);
sprintf(jac,"%s_s_LBFGS_rho_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_rho_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]);
remove(jac);
}
if(LBFGS_NPAR>2){
fclose(FP_JAC3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_s_LBFGS_vp_it%d_w%d.bin",JACOBIAN,iter+1,w);
if (MYID==0) mergemod(jac,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(jac,"%s_bfgsmod_rho_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]);
sprintf(jac,"%s_s_LBFGS_vp_it%d_w%d.bin.%i.%i",JACOBIAN,iter+1,w,POS[1],POS[2]);
remove(jac);
}
}
......
src/denise.c
View file @ e85740fe
......@@ -108,10 +108,10 @@ int main(int argc, char **argv){
int QUELLART_OLD;
/* Variables for L-BFGS */
int LBFGS=0, bfgsnum=0,bfgspar=3;
int LBFGS=0,LBFGS_NPAR=3;
int LBFGS_iter_start=1;
float LBFGS_L2_temp;
float **bfgsmod1,**bfgsmod2,**bfgsmod3,**bfgsgrad1, **bfgsgrad2, **bfgsgrad3, *bfgsscale1,*bfgsscale2,*bfgsscale3;
float **s_LBFGS,**bfgsmod2,**bfgsmod3,**y_LBFGS, **bfgsgrad2, **bfgsgrad3, *rho_LBFGS,*bfgsscale2,*bfgsscale3;
int l=0;
int w=0;
int m=0;
......@@ -391,22 +391,22 @@ int main(int argc, char **argv){
}
if(GRAD_METHOD==2) {
bfgsnum=N_LBFGS;
/* Allocate memory for L-BFGS */
if(WAVETYPE==2) bfgspar=2;
if(WAVETYPE==2) LBFGS_NPAR=2;
bfgsmod1=fmatrix(1,bfgsnum,1,bfgspar*NX*NY);
s_LBFGS=fmatrix(1,N_LBFGS,1,LBFGS_NPAR*NX*NY);
bfgsgrad1=fmatrix(1,bfgsnum,1,bfgspar*NX*NY);
y_LBFGS=fmatrix(1,N_LBFGS,1,LBFGS_NPAR*NX*NY);
bfgsscale1=vector(1,bfgsnum);
rho_LBFGS=vector(1,N_LBFGS);
for(l=1;l<=bfgsnum;l++){
for(m=1;m<=bfgspar*NX*NY;m++){
bfgsmod1[l][m]=0.0;
bfgsgrad1[l][m]=0.0;
for(l=1;l<=N_LBFGS;l++){
for(m=1;m<=LBFGS_NPAR*NX*NY;m++){
s_LBFGS[l][m]=0.0;
y_LBFGS[l][m]=0.0;
}
bfgsscale1[l]=0.0;
rho_LBFGS[l]=0.0;
}
}
......@@ -969,7 +969,7 @@ int main(int argc, char **argv){
/* restart L-BFGS */
if(iter==LBFGS_iter_start) {
lbfgs_reset(iter,bfgsnum,bfgspar,bfgsmod1,bfgsgrad1,bfgsscale1);
lbfgs_reset(iter,N_LBFGS,LBFGS_NPAR,s_LBFGS,y_LBFGS,rho_LBFGS);
/* set values */
FWI_run=1;
......@@ -2913,7 +2913,7 @@ int main(int argc, char **argv){
/*---------------------*/
/* L-BFGS */
/*---------------------*/
lbfgs(waveconv_u, waveconv_rho, waveconv,Vs_avg,rho_avg,Vp_avg, bfgsscale1, bfgsmod1, bfgsgrad1, bfgsnum,bfgspar, iter,&LBFGS_iter_start);
lbfgs(waveconv_u, waveconv_rho, waveconv,Vs_avg,rho_avg,Vp_avg, rho_LBFGS, s_LBFGS, y_LBFGS, N_LBFGS,LBFGS_NPAR, iter,&LBFGS_iter_start);
}
......@@ -2961,14 +2961,14 @@ int main(int argc, char **argv){
if(iter>2 && WOLFE_TRY_OLD_STEPLENGTH) alpha_SL=alpha_SL_old;
/* Calculate update */
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,INVMAT,alpha_SL,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,bfgsmod1,bfgsnum,bfgspar,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,INVMAT,alpha_SL,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
FWI_run=1;
}
/* Check if current step length satisfy wolfe condition, if not call linesearch for new step length */
if(countstep>0) {
wolfe_status=check_wolfe(alpha_SL, L2_SL_old, L2_SL_new, waveconv_u_old, waveconv_u, waveconv_u_up,waveconv_rho_old, waveconv_rho, waveconv_rho_up,waveconv_old, waveconv, waveconv_up, c1_SL, c2_SL,bfgspar);
wolfe_status=check_wolfe(alpha_SL, L2_SL_old, L2_SL_new, waveconv_u_old, waveconv_u, waveconv_u_up,waveconv_rho_old, waveconv_rho, waveconv_rho_up,waveconv_old, waveconv, waveconv_up, c1_SL, c2_SL,LBFGS_NPAR);
if(wolfe_status==0) {
/* Current step length satisfy wolfe condition, abort step length search */
......@@ -2989,7 +2989,7 @@ int main(int argc, char **argv){
/* Current step length do not satisfy wolfe condition, try new step length */
wolfe_linesearch(wolfe_status, &alpha_SL_min, &alpha_SL_max, &alpha_SL);
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,INVMAT,alpha_SL,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,bfgsmod1,bfgsnum,bfgspar,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,INVMAT,alpha_SL,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
fprintf(FP,"; New steplength=%1.3f",alpha_SL);
FWI_run=1;
}
......@@ -3023,7 +3023,7 @@ int main(int argc, char **argv){
}
/* do the final model update */
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,bfgsmod1,bfgsnum,bfgspar,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
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);
/* write L2 log file */
if (TIME_FILT==0){
......@@ -3102,7 +3102,7 @@ int main(int argc, char **argv){
for (itest=itests;itest<=iteste;itest++){
/* calculate change in the material parameters */
calc_mat_change_test(waveconv,waveconv_rho,waveconv_u,prho,prhonp1,ppi,ppinp1,pu,punp1,iter,1,INVMAT,eps_scale,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,bfgsmod1,bfgsnum,bfgspar,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
calc_mat_change_test(waveconv,waveconv_rho,waveconv_u,prho,prhonp1,ppi,ppinp1,pu,punp1,iter,1,INVMAT,eps_scale,1,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
/* For the calculation of the material parameters beteween gridpoints
the have to be averaged. For this, values lying at 0 and NX+1,
......@@ -3660,7 +3660,7 @@ int main(int argc, char **argv){
/* -----------------------------------------------------------------------*/
/* ----------- Do the actual update to the material parameters -----------*/
/* -----------------------------------------------------------------------*/
calc_mat_change_test(waveconv,waveconv_rho,waveconv_u,prho,prhonp1,ppi,ppinp1,pu,punp1,iter,1,INVMAT,eps_scale,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,bfgsmod1,bfgsnum,bfgspar,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
calc_mat_change_test(waveconv,waveconv_rho,waveconv_u,prho,prhonp1,ppi,ppinp1,pu,punp1,iter,1,INVMAT,eps_scale,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
fprintf(FP,"=================================================================================================\n");
} /* end of if(INVMAT!=4) */
......
src/read_par_json.c
View file @ e85740fe
......@@ -708,6 +708,9 @@ void read_par_json(FILE *fp, char *fileinp){
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);
......@@ -749,8 +752,6 @@ void read_par_json(FILE *fp, char *fileinp){
}
/* Trace killing STF */
if (get_int_from_objectlist("TRKILL_STF",number_readobjects,&TRKILL_STF,varname_list, value_list)){
TRKILL_STF=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