Commit c7b43d05 authored by Florian Wittkamp's avatar Florian Wittkamp

Maintenance of the Code

(1.)
Renamed function err() to declare_error().
This was a problem, as many standard libraries also has functions err().

(2)
I fixed a lot of compiler warnings like unused variables and uninitialised variables. However, here is a lot of work to do.
parent 03134342
...@@ -46,7 +46,7 @@ void model_acoustic(float ** rho, float ** pi){ ...@@ -46,7 +46,7 @@ void model_acoustic(float ** rho, float ** pi){
/*-----------------------------------------------------------------------*/ /*-----------------------------------------------------------------------*/
y=h/DH; y=h/DH;
if(y==NYG) err(" \n y is equal NYG !! see src/model_grad.c \n "); if(y==NYG) declare_error(" \n y is equal NYG !! see src/model_grad.c \n ");
grad1=(vp2-vp1)/y; grad1=(vp2-vp1)/y;
grad3=(rho2-rho1)/y; grad3=(rho2-rho1)/y;
......
...@@ -46,7 +46,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){ ...@@ -46,7 +46,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){
/*-----------------------------------------------------------------------*/ /*-----------------------------------------------------------------------*/
y=h/DH; y=h/DH;
if(y==NYG) err(" \n y is equal NYG !! see src/model_grad.c \n "); if(y==NYG) declare_error(" \n y is equal NYG !! see src/model_grad.c \n ");
grad1=(vp2-vp1)/y; grad1=(vp2-vp1)/y;
grad2=(vs2-vs1)/y; grad2=(vs2-vs1)/y;
grad3=(rho2-rho1)/y; grad3=(rho2-rho1)/y;
......
...@@ -52,7 +52,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float ** ...@@ -52,7 +52,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float **
y=h/DH; y=h/DH;
if(y==NYG) err(" \n y is equal NYG !! see src/model_grad.c \n "); if(y==NYG) declare_error(" \n y is equal NYG !! see src/model_grad.c \n ");
grad1=(vp2-vp1)/y; grad1=(vp2-vp1)/y;
grad2=(vs2-vs1)/y; grad2=(vs2-vs1)/y;
grad3=(rho2-rho1)/y; grad3=(rho2-rho1)/y;
......
...@@ -52,7 +52,7 @@ void model_viscac(float ** rho, float ** pi, float ** taup, float * eta){ ...@@ -52,7 +52,7 @@ void model_viscac(float ** rho, float ** pi, float ** taup, float * eta){
y=h/DH; y=h/DH;
if(y==NYG) err(" \n y is equal NYG !! see src/model_grad.c \n "); if(y==NYG) declare_error(" \n y is equal NYG !! see src/model_grad.c \n ");
grad1=(vp2-vp1)/y; grad1=(vp2-vp1)/y;
grad3=(rho2-rho1)/y; grad3=(rho2-rho1)/y;
......
...@@ -38,7 +38,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){ ...@@ -38,7 +38,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){
/* local variables */ /* local variables */
float rhov, muv, piv, vp, vs; float rhov, vp, vs;
float vp0, vs0, rho0, gvp, gvs, grho; float vp0, vs0, rho0, gvp, gvs, grho;
int i, j, ii, jj; int i, j, ii, jj;
FILE *fp_vs, *fp_vp, *fp_rho; FILE *fp_vs, *fp_vp, *fp_rho;
...@@ -66,18 +66,18 @@ void model_elastic(float ** rho, float ** pi, float ** u){ ...@@ -66,18 +66,18 @@ void model_elastic(float ** rho, float ** pi, float ** u){
fprintf(FP,"\t Vp:\n\t %s.vp\n\n",MFILE); fprintf(FP,"\t Vp:\n\t %s.vp\n\n",MFILE);
sprintf(filename,"%s.vp",MFILE); sprintf(filename,"%s.vp",MFILE);
fp_vp=fopen(filename,"r"); fp_vp=fopen(filename,"r");
if (fp_vp==NULL) err(" Could not open model file for Vp ! "); if (fp_vp==NULL) declare_error(" Could not open model file for Vp ! ");
fprintf(FP,"\t Vs:\n\t %s.vs\n\n",MFILE); fprintf(FP,"\t Vs:\n\t %s.vs\n\n",MFILE);
sprintf(filename,"%s.vs",MFILE); sprintf(filename,"%s.vs",MFILE);
fp_vs=fopen(filename,"r"); fp_vs=fopen(filename,"r");
if (fp_vs==NULL) err(" Could not open model file for Vs ! "); if (fp_vs==NULL) declare_error(" Could not open model file for Vs ! ");
fprintf(FP,"\t Density:\n\t %s.rho\n\n",MFILE); fprintf(FP,"\t Density:\n\t %s.rho\n\n",MFILE);
sprintf(filename,"%s.rho",MFILE); sprintf(filename,"%s.rho",MFILE);
fp_rho=fopen(filename,"r"); fp_rho=fopen(filename,"r");
if (fp_rho==NULL) err(" Could not open model file for densities ! "); if (fp_rho==NULL) declare_error(" Could not open model file for densities ! ");
} }
/* read density and Lame parameters */ /* read density and Lame parameters */
...@@ -86,18 +86,18 @@ void model_elastic(float ** rho, float ** pi, float ** u){ ...@@ -86,18 +86,18 @@ void model_elastic(float ** rho, float ** pi, float ** u){
fprintf(FP,"\t Lame parameter lambda:\n\t %s.lam\n\n",MFILE); fprintf(FP,"\t Lame parameter lambda:\n\t %s.lam\n\n",MFILE);
sprintf(filename,"%s.lam",MFILE); sprintf(filename,"%s.lam",MFILE);
fp_vp=fopen(filename,"r"); fp_vp=fopen(filename,"r");
if (fp_vp==NULL) err(" Could not open model file for Lame parameter lambda ! "); if (fp_vp==NULL) declare_error(" Could not open model file for Lame parameter lambda ! ");
fprintf(FP,"\t Lame parameter mu:\n\t %s.vs\n\n",MFILE); fprintf(FP,"\t Lame parameter mu:\n\t %s.vs\n\n",MFILE);
sprintf(filename,"%s.mu",MFILE); sprintf(filename,"%s.mu",MFILE);
fp_vs=fopen(filename,"r"); fp_vs=fopen(filename,"r");
if (fp_vs==NULL) err(" Could not open model file for Lame parameter mu ! "); if (fp_vs==NULL) declare_error(" Could not open model file for Lame parameter mu ! ");
fprintf(FP,"\t Density:\n\t %s.rho\n\n",MFILE); fprintf(FP,"\t Density:\n\t %s.rho\n\n",MFILE);
sprintf(filename,"%s.rho",MFILE); sprintf(filename,"%s.rho",MFILE);
fp_rho=fopen(filename,"r"); fp_rho=fopen(filename,"r");
if (fp_rho==NULL) err(" Could not open model file for densities ! "); if (fp_rho==NULL) declare_error(" Could not open model file for densities ! ");
} }
......
...@@ -67,7 +67,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float ** ...@@ -67,7 +67,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float **
/*read FL nodes from File*/ /*read FL nodes from File*/
flfile=fopen("model/final.mod.flnodes.Q20","r"); flfile=fopen("model/final.mod.flnodes.Q20","r");
if (flfile==NULL) err(" FL-file could not be opened !"); if (flfile==NULL) declare_error(" FL-file could not be opened !");
......
...@@ -33,7 +33,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float ** ...@@ -33,7 +33,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float **
extern char INV_MODELFILE[STRING_SIZE]; extern char INV_MODELFILE[STRING_SIZE];
extern float DH, *FL, TAU, DT; extern float DH, *FL, TAU, DT;
/* local variables */ /* local variables */
float vp, vs, rhov, grad, y1, y2, ts, tp, muv, piv, *pts; float vp, vs, rhov, grad, y1, y2, ts, tp, *pts;
int i, j, ii, jj, l; int i, j, ii, jj, l;
char modfile[STRING_SIZE]; char modfile[STRING_SIZE];
...@@ -53,7 +53,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float ** ...@@ -53,7 +53,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float **
y1=h/DH; y1=h/DH;
y2=layer/DH; y2=layer/DH;
if(y1==NYG) err(" \n y is equal NYG !! see src/model_grad.c \n "); if(y1==NYG) declare_error(" \n y is equal NYG !! see src/model_grad.c \n ");
grad=(vs2-vs1)/y1; grad=(vs2-vs1)/y1;
......
...@@ -57,7 +57,7 @@ void model_acoustic(float **rho, float **pi){ ...@@ -57,7 +57,7 @@ void model_acoustic(float **rho, float **pi){
if(SWS_TAPER_FILE) taper=matrix(-nd+1,NY+nd,-nd+1,NX+nd); if(SWS_TAPER_FILE) taper=matrix(-nd+1,NY+nd,-nd+1,NX+nd);
flfile=fopen("model_true/flnodes.toy_example_ac.start","r"); flfile=fopen("model_true/flnodes.toy_example_ac.start","r");
if (flfile==NULL) err(" FL-file could not be opened !"); if (flfile==NULL) declare_error(" FL-file could not be opened !");
/* Read parameters */ /* Read parameters */
for (l=1;l<=nodes;l++){ for (l=1;l<=nodes;l++){
...@@ -187,7 +187,7 @@ void model_acoustic(float **rho, float **pi){ ...@@ -187,7 +187,7 @@ void model_acoustic(float **rho, float **pi){
flrho=vector(1,nodes); flrho=vector(1,nodes);
flvp=vector(1,nodes); flvp=vector(1,nodes);
flfile=fopen("model_true/flnodes.toy_example_ac","r"); flfile=fopen("model_true/flnodes.toy_example_ac","r");
if (flfile==NULL) err(" FL-file could not be opened !"); if (flfile==NULL) declare_error(" FL-file could not be opened !");
/* Read parameters */ /* Read parameters */
for (l=1;l<=nodes;l++){ for (l=1;l<=nodes;l++){
......
...@@ -63,7 +63,7 @@ void model_acoustic(float **rho, float **pi){ ...@@ -63,7 +63,7 @@ void model_acoustic(float **rho, float **pi){
/*read FL nodes from File*/ /*read FL nodes from File*/
flfile=fopen("model_true/flnodes.toy_example_ac","r"); flfile=fopen("model_true/flnodes.toy_example_ac","r");
if (flfile==NULL) err(" FL-file could not be opened !"); if (flfile==NULL) declare_error(" FL-file could not be opened !");
for (l=1;l<=nodes;l++){ for (l=1;l<=nodes;l++){
......
...@@ -38,7 +38,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){ ...@@ -38,7 +38,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){
extern char INV_MODELFILE[STRING_SIZE]; extern char INV_MODELFILE[STRING_SIZE];
/* local variables */ /* local variables */
float muv, piv, vp, vs, rhov, ts, tp, *pts; float muv, piv, vp, vs, rhov, ts, tp;
int i, j, ii, jj, l; int i, j, ii, jj, l;
char modfile[STRING_SIZE]; char modfile[STRING_SIZE];
...@@ -62,7 +62,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){ ...@@ -62,7 +62,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){
flrho=vector(1,nodes); flrho=vector(1,nodes);
flvp=vector(1,nodes); flvp=vector(1,nodes);
flfile=fopen("model_true/flnodes.toy_example.start","r"); flfile=fopen("model_true/flnodes.toy_example.start","r");
if (flfile==NULL) err(" FL-file could not be opened !"); if (flfile==NULL) declare_error(" FL-file could not be opened !");
...@@ -154,7 +154,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){ ...@@ -154,7 +154,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){
flrho=vector(1,nodes); flrho=vector(1,nodes);
flvp=vector(1,nodes); flvp=vector(1,nodes);
flfile=fopen("model_true/flnodes.toy_example","r"); flfile=fopen("model_true/flnodes.toy_example","r");
if (flfile==NULL) err(" FL-file could not be opened !"); if (flfile==NULL) declare_error(" FL-file could not be opened !");
for (l=1;l<=nodes;l++){ for (l=1;l<=nodes;l++){
fgets(cline,255,flfile); fgets(cline,255,flfile);
......
...@@ -55,7 +55,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){ ...@@ -55,7 +55,7 @@ void model_elastic(float ** rho, float ** pi, float ** u){
/*read FL nodes from File*/ /*read FL nodes from File*/
flfile=fopen("model_true/flnodes.toy_example","r"); flfile=fopen("model_true/flnodes.toy_example","r");
if (flfile==NULL) err(" FL-file could not be opened !"); if (flfile==NULL) declare_error(" FL-file could not be opened !");
......
...@@ -67,7 +67,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float ** ...@@ -67,7 +67,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float **
flrho=vector(1,nodes); flrho=vector(1,nodes);
flvp=vector(1,nodes); flvp=vector(1,nodes);
flfile=fopen("model_true/flnodes.toy_example.start","r"); flfile=fopen("model_true/flnodes.toy_example.start","r");
if (flfile==NULL) err(" FL-file could not be opened !"); if (flfile==NULL) declare_error(" FL-file could not be opened !");
...@@ -165,7 +165,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float ** ...@@ -165,7 +165,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float **
flrho=vector(1,nodes); flrho=vector(1,nodes);
flvp=vector(1,nodes); flvp=vector(1,nodes);
flfile=fopen("model_true/flnodes.toy_example","r"); flfile=fopen("model_true/flnodes.toy_example","r");
if (flfile==NULL) err(" FL-file could not be opened !"); if (flfile==NULL) declare_error(" FL-file could not be opened !");
for (l=1;l<=nodes;l++){ for (l=1;l<=nodes;l++){
fgets(cline,255,flfile); fgets(cline,255,flfile);
......
...@@ -65,7 +65,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float ** ...@@ -65,7 +65,7 @@ void model(float ** rho, float ** pi, float ** u, float ** taus, float **
/*read FL nodes from File*/ /*read FL nodes from File*/
flfile=fopen("model_true/flnodes.toy_example","r"); flfile=fopen("model_true/flnodes.toy_example","r");
if (flfile==NULL) err(" FL-file could not be opened !"); if (flfile==NULL) declare_error(" FL-file could not be opened !");
......
...@@ -32,24 +32,24 @@ ...@@ -32,24 +32,24 @@
int main(int argc, char **argv){ int main(int argc, char **argv){
/* variables in main */ /* variables in main */
int ns, nseismograms=0, nt, nd, fdo3, j, i, ii, jj, shotid, recid, k, nc, iter, h, infoout, SHOTINC, TIMEWIN, test_eps, lq, iq, jq, hin, hin1, s=0; int ns, nseismograms=0, nt, nd, fdo3, j, i, iter, h, infoout, SHOTINC, hin, hin1, s=0;
int NTDTINV, nxny, nxnyi, imat, imat1, imat2, IDXI, IDYI, hi, NTST, NTSTI, partest; int NTDTINV, nxny, nxnyi, imat, imat1, imat2, IDXI, IDYI, hi, NTST, NTSTI;
int lsnap, nsnap=0, lsamp=0, buffsize, invtime, invtimer, sws, swstestshot, snapseis, snapseis1, PML; int lsnap, nsnap=0, lsamp=0, buffsize, swstestshot, snapseis, snapseis1;
int ntr=0, ntr_loc=0, ntr_glob=0, nsrc=0, nsrc_loc=0, nsrc_glob=0, ishot, irec, nshots=0, nshots1, Lcount, itest, Lcountsum, itestshot; int ntr=0, ntr_loc=0, ntr_glob=0, nsrc=0, nsrc_loc=0, nsrc_glob=0, ishot, irec, nshots=0, nshots1, Lcount, itest, itestshot;
float pum, ppim, ppim1, ppim2, thetaf, thetab, e33, e33b, e11, e11b, muss, lamss; float muss, lamss;
float memdyn, memmodel, memseismograms, membuffer, memtotal, dngn, fphi, sum, avggrad, beta, betan, betaz, betaLog, betaVp, betaVs, betarho, eps_scale, L2old; float memdyn, memmodel, memseismograms, membuffer, memtotal, eps_scale;
float fac1, fac2, wavefor, waverecipro, dump, dump1, epsilon, gradsign, mun, eps1, gradplastiter, gradglastiter, gradclastiter, betar, sig_max, sig_max1; float fac1, fac2;
float signL1, RMS, opteps_vp, opteps_vs, opteps_rho, Vs, Vp, Vp_avg, C_vp, Vs_avg, C_vs, Cd, rho_avg, C_rho, Vs_sum, Vp_sum, rho_sum, Zp, Zs; float opteps_vp, opteps_vs, opteps_rho, Vp_avg, C_vp, Vs_avg, C_vs, rho_avg, C_rho;
float freqshift, dfreqshift, memfwt, memfwt1, memfwtdata; float memfwt, memfwt1, memfwtdata;
char *buff_addr, ext[10], *fileinp; char *buff_addr, ext[10], *fileinp;
char wave_forward[225], wave_recipro[225], wave_conv[225], jac[225], jac2[225], jacsum[225], dwavelet[225], vyf[STRING_SIZE]; char jac[225];
double time1, time2, time3, time4, time5, time6, time7, time8, double time1, time2, time3, time4, time5, time6, time7, time8,
time_av_v_update=0.0, time_av_s_update=0.0, time_av_v_exchange=0.0, time_av_v_update=0.0, time_av_s_update=0.0, time_av_v_exchange=0.0,
time_av_s_exchange=0.0, time_av_timestep=0.0; time_av_s_exchange=0.0, time_av_timestep=0.0;
float L2, L2sum, L2_all_shots, L2sum_all_shots, *L2t, alphanomsum, alphanom, alphadenomsum, alphadenom, scaleamp ,sdummy, lamr; float L2, L2sum, L2_all_shots, L2sum_all_shots, *L2t, alphanom, alphadenom;
int sum_killed_traces=0, sum_killed_traces_testshots=0, killed_traces=0, killed_traces_testshots=0; int sum_killed_traces=0, sum_killed_traces_testshots=0, killed_traces=0, killed_traces_testshots=0;
int *ptr_killed_traces=&killed_traces, *ptr_killed_traces_testshots=&killed_traces_testshots; int *ptr_killed_traces=&killed_traces, *ptr_killed_traces_testshots=&killed_traces_testshots;
...@@ -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, ** uz, ** uxy, ** uyx, ** Vp0, ** uttx, ** utty, ** Vs0, ** Rho0; float ** psxx, ** psxy, ** psyy, ** psxz, ** psyz, **psp, ** ux, ** uy, ** uxy, ** uyx, ** 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;
...@@ -74,7 +74,7 @@ int main(int argc, char **argv){ ...@@ -74,7 +74,7 @@ int main(int argc, char **argv){
** sectionvydiff=NULL, ** sectionpn=NULL, ** sectionread=NULL, ** sectionvy_conv=NULL, ** sectionvy_obs=NULL, ** sectionvx_conv=NULL,** sectionvx_obs=NULL, ** sectionvz_conv=NULL,** sectionvz_obs=NULL, ** sectionvydiff=NULL, ** sectionpn=NULL, ** sectionread=NULL, ** sectionvy_conv=NULL, ** sectionvy_obs=NULL, ** sectionvx_conv=NULL,** sectionvx_obs=NULL, ** sectionvz_conv=NULL,** sectionvz_obs=NULL,
** sectionp_conv=NULL,** sectionp_obs=NULL, * source_time_function=NULL; ** sectionp_conv=NULL,** sectionp_obs=NULL, * source_time_function=NULL;
float ** absorb_coeff, ** taper_coeff, * epst1, * epst2, * epst3, * picked_times; float ** absorb_coeff, ** taper_coeff, * epst1, * epst2, * epst3, * picked_times;
float ** srcpos=NULL, **srcpos_loc=NULL, ** srcpos1=NULL, **srcpos_loc_back=NULL, ** signals=NULL,** signals_SH=NULL, ** signals_rec=NULL, *hc=NULL; float ** srcpos=NULL, **srcpos_loc=NULL, ** srcpos1=NULL, **srcpos_loc_back=NULL, ** signals=NULL,** signals_SH=NULL, *hc=NULL;
int ** recpos=NULL, ** recpos_loc=NULL; int ** recpos=NULL, ** recpos_loc=NULL;
/*int ** tracekill=NULL, TRKILL, DTRKILL;*/ /*int ** tracekill=NULL, TRKILL, DTRKILL;*/
int * DTINV_help; int * DTINV_help;
...@@ -100,12 +100,10 @@ int main(int argc, char **argv){ ...@@ -100,12 +100,10 @@ int main(int argc, char **argv){
int SOURCE_SHAPE_OLD; int SOURCE_SHAPE_OLD;
/* Variables for L-BFGS */ /* Variables for L-BFGS */
int LBFGS=0,LBFGS_NPAR=3; int LBFGS_NPAR=3;
int LBFGS_iter_start=1; int LBFGS_iter_start=1;
float LBFGS_L2_temp;
float **s_LBFGS,**y_LBFGS, *rho_LBFGS; float **s_LBFGS,**y_LBFGS, *rho_LBFGS;
int l=0; int l=0;
int w=0;
int m=0; int m=0;
/* Check wolfe */ /* Check wolfe */
...@@ -136,9 +134,6 @@ int main(int argc, char **argv){ ...@@ -136,9 +134,6 @@ int main(int argc, char **argv){
int * recswitch=NULL; int * recswitch=NULL;
float ** fulldata=NULL, ** fulldata_vx=NULL, ** fulldata_vy=NULL, ** fulldata_vz=NULL, ** fulldata_p=NULL, ** fulldata_curl=NULL, ** fulldata_div=NULL; float ** fulldata=NULL, ** fulldata_vx=NULL, ** fulldata_vy=NULL, ** fulldata_vz=NULL, ** fulldata_p=NULL, ** fulldata_curl=NULL, ** fulldata_div=NULL;
/* different modelling types */
int mod_type=0;
/*vector for abort criterion*/ /*vector for abort criterion*/
float * L2_hist=NULL; float * L2_hist=NULL;
...@@ -158,11 +153,9 @@ int main(int argc, char **argv){ ...@@ -158,11 +153,9 @@ int main(int argc, char **argv){
float *FC_EXT=NULL; float *FC_EXT=NULL;
int nfrq=0; int nfrq=0;
int FREQ_NR=1; int FREQ_NR=1;
/* declaration of variables for trace killing */
int ** kill_tmp;
FILE *ftracekill;
FILE *fprec, *FP2, *FP3, *FP4, *FP5, *FPL2, *FP6, *FP7; FILE *fprec, *FPL2;
/* General parameters */ /* General parameters */
int nt_out; int nt_out;
...@@ -193,7 +186,7 @@ int main(int argc, char **argv){ ...@@ -193,7 +186,7 @@ int main(int argc, char **argv){
printf("\n==================================================================\n"); printf("\n==================================================================\n");
printf(" Cannot open IFOS input file %s \n",fileinp); printf(" Cannot open IFOS input file %s \n",fileinp);
printf("\n==================================================================\n\n"); printf("\n==================================================================\n\n");
err(" --- "); declare_error(" --- ");
} }
} }
...@@ -346,7 +339,7 @@ int main(int argc, char **argv){ ...@@ -346,7 +339,7 @@ int main(int argc, char **argv){
/* allocate buffer for buffering messages */ /* allocate buffer for buffering messages */
buff_addr=malloc(buffsize); buff_addr=malloc(buffsize);
if (!buff_addr) err("allocation failure for buffer for MPI_Bsend !"); if (!buff_addr) declare_error("allocation failure for buffer for MPI_Bsend !");
MPI_Buffer_attach(buff_addr,buffsize); MPI_Buffer_attach(buff_addr,buffsize);
/* allocation for request and status arrays */ /* allocation for request and status arrays */
...@@ -1196,7 +1189,7 @@ int main(int argc, char **argv){ ...@@ -1196,7 +1189,7 @@ int main(int argc, char **argv){
srcpos_loc = splitsrc(srcpos,&nsrc_loc, nsrc); srcpos_loc = splitsrc(srcpos,&nsrc_loc, nsrc);
} }
if((SOURCE_SHAPE==7)||(SOURCE_SHAPE==3))err("SOURCE_SHAPE==7 or SOURCE_SHAPE==3 isn't possible with INV_STF==1"); if((SOURCE_SHAPE==7)||(SOURCE_SHAPE==3))declare_error("SOURCE_SHAPE==7 or SOURCE_SHAPE==3 isn't possible with INV_STF==1");
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
...@@ -1256,7 +1249,7 @@ int main(int argc, char **argv){ ...@@ -1256,7 +1249,7 @@ int main(int argc, char **argv){
if (WAVETYPE==1 || WAVETYPE==3) { if (WAVETYPE==1 || WAVETYPE==3) {
if (isnan(pvy[NY/2][NX/2])) { if (isnan(pvy[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]); fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !"); declare_error(" Simulation is unstable !");
} }
} }
...@@ -1264,7 +1257,7 @@ int main(int argc, char **argv){ ...@@ -1264,7 +1257,7 @@ int main(int argc, char **argv){
if (WAVETYPE==2 || WAVETYPE==3) { if (WAVETYPE==2 || WAVETYPE==3) {
if (isnan(pvz[NY/2][NX/2])) { if (isnan(pvz[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]); fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !"); declare_error(" Simulation is unstable !");
} }
} }
...@@ -1685,7 +1678,7 @@ int main(int argc, char **argv){ ...@@ -1685,7 +1678,7 @@ int main(int argc, char **argv){
if (WAVETYPE==1 || WAVETYPE==3) { if (WAVETYPE==1 || WAVETYPE==3) {
if (isnan(pvy[NY/2][NX/2])) { if (isnan(pvy[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]); fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !"); declare_error(" Simulation is unstable !");
} }
} }
...@@ -1693,7 +1686,7 @@ int main(int argc, char **argv){ ...@@ -1693,7 +1686,7 @@ int main(int argc, char **argv){
if (WAVETYPE==2 || WAVETYPE==3) { if (WAVETYPE==2 || WAVETYPE==3) {
if (isnan(pvz[NY/2][NX/2])) { if (isnan(pvz[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]); fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !"); declare_error(" Simulation is unstable !");
} }
} }
...@@ -2187,14 +2180,14 @@ int main(int argc, char **argv){ ...@@ -2187,14 +2180,14 @@ int main(int argc, char **argv){
if (WAVETYPE==1 || WAVETYPE==3) { if (WAVETYPE==1 || WAVETYPE==3) {
if (isnan(pvy[NY/2][NX/2])) { if (isnan(pvy[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]); fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !"); declare_error(" Simulation is unstable !");
} }
} }
/* Check if simulation is still stable SH */ /* Check if simulation is still stable SH */
if (WAVETYPE==2 || WAVETYPE==3) { if (WAVETYPE==2 || WAVETYPE==3) {
if (isnan(pvz[NY/2][NX/2])) { if (isnan(pvz[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]); fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !"); declare_error(" Simulation is unstable !");
} }
} }
...@@ -3332,14 +3325,14 @@ int main(int argc, char **argv){ ...@@ -3332,14 +3325,14 @@ int main(int argc, char **argv){
if (WAVETYPE==1 || WAVETYPE==3) { if (WAVETYPE==1 || WAVETYPE==3) {
if (isnan(pvy[NY/2][NX/2])) { if (isnan(pvy[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]); fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !"); declare_error(" Simulation is unstable !");
} }
} }
/* Check if simulation is still stable SH */ /* Check if simulation is still stable SH */
if (WAVETYPE==2 || WAVETYPE==3) { if (WAVETYPE==2 || WAVETYPE==3) {
if (isnan(pvz[NY/2][NX/2])) { if (isnan(pvz[NY/2][NX/2])) {
fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]); fprintf(FP,"\n Time step: %d; pvy: %f \n",nt,pvy[NY/2][NX/2]);
err(" Simulation is unstable !"); declare_error(" Simulation is unstable !");
} }
} }
......
...@@ -36,14 +36,11 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float ...@@ -36,14 +36,11 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
extern int ACOUSTIC; extern int ACOUSTIC;
/* local */ /* local */
int m=0,v=0,w=0; int m=0,w=0;
int i,j,k,l; int i,j,l;
float beta_LBFGS=0.0;
float dum1=0.0, dum2=0.0, buf1=0.0, buf2=0.0;
float *q_LBFGS,*alpha_LBFGS,*r_LBFGS; float *q_LBFGS,*alpha_LBFGS,*r_LBFGS;
float h0;
char jac[225]; char jac[225];
FILE *FP_JAC; FILE *FP_JAC = NULL;
/*---------------------*/ /*---------------------*/
...@@ -90,9 +87,6 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float ...@@ -90,9 +87,6 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
q_LBFGS = vector(1,NPAR_LBFGS*NX*NY); q_LBFGS = vector(1,NPAR_LBFGS*NX*NY);
r_LBFGS = vector(1,NPAR_LBFGS*NX*NY); r_LBFGS = vector(1,NPAR_LBFGS*NX*NY);
m=iteration-N_LBFGS;
if(m<1) m=1;
w=(iteration-1)%N_LBFGS; w=(iteration-1)%N_LBFGS;
if(w==0) w=N_LBFGS; if(w==0) w=N_LBFGS;
...@@ -217,7 +211,6 @@ void lbfgs_core(int iteration, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS, flo ...@@ -217,7 +211,6 @@ void lbfgs_core(int iteration, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS, flo
float dum1=0.0, dum2=0.0, buf1=0.0, buf2=0.0; float dum1=0.0, dum2=0.0, buf1=0.0, buf2=0.0;
float h0; float h0;
int m=0,v=0,w=0,l=0; int m=0,v=0,w=0,l=0;
int VERBOSE_local=1;
m=iteration-N_LBFGS; m=iteration-N_LBFGS;
if(m<1) m=1; if(m<1) m=1;
...@@ -253,7 +246,6 @@ void lbfgs_core(int iteration, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS, flo ...@@ -253,7 +246,6 @@ void lbfgs_core(int iteration, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS, flo
/* L-BFGS loop 1 */ /* L-BFGS loop 1 */
/*----------------------------------*/ /*----------------------------------*/
l=0;
for(v=iteration-1; v>=m;v--){ for(v=iteration-1; v>=m;v--){
w=v%N_LBFGS; w=v%N_LBFGS;
......
...@@ -39,7 +39,7 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ...@@ -39,7 +39,7 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int
extern FILE *FP; extern FILE *FP;
int use_conjugate_1=1; int use_conjugate_1=1;
int use_conjugate_2=1; int use_conjugate_2=1;
FILE *FP3, *FP4, *FP6, *FP5; FILE *FP3, *FP4, *FP6, *FP5 = NULL;
...@@ -47,7 +47,7 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ...@@ -47,7 +47,7 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int
/* ===================================================== GRADIENT ZP ================================================================================== */ /* ===================================================== GRADIENT ZP ================================================================================== */
/* =========================================================================