Skip to content

IFOS2D
IFOS2D
Commit 0904becf authored by Florian Wittkamp's avatar Florian Wittkamp
Browse files
Options

Offset based TK now additionally to TK File 

If TRKILL_STF_OFFSET or TRKILL_OFFSET is set to 2, the offset based TraceKill will be applied in addition to the TraceKill file.
Therefore, it is possible to exclude with TraceKill file fault receiveres and use an offset based tracekill with the build in feature at the same time.
parent b57d6b66
Hide whitespace changes
Inline Side-by-side
Showing with 154 additions and 43 deletions
src/LBFGS.c
View file @ 0904becf
......@@ -34,9 +34,10 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
extern char JACOBIAN[STRING_SIZE];
extern int WAVETYPE;
extern int ACOUSTIC;
extern float LBFGS_SCALE_GRADIENTS;
/* local */
int m=0,w=0;
int w=0;
int i,j,l;
float *q_LBFGS,*alpha_LBFGS,*r_LBFGS;
char jac[225];
......@@ -60,6 +61,26 @@ void lbfgs(float **grad_vs, float **grad_rho, float **grad_vp,float Vs_avg,float
write_matrix_disk(grad_vp, jac);
}
/*---------------------*/
/* Experimental */
/*---------------------*/
/* Scale the gradients with a constant factor */
/* This is to avoid numerical instabilities if the gradients are to small (absolute value) */
/* Do not use this feature, unless you have to */
if(LBFGS_SCALE_GRADIENTS!=1) {
if(MYID==0) printf("\n\n Scaling the gradients to ensure L-BFGS stability");
if(MYID==0) printf("\n Scaling with %f",LBFGS_SCALE_GRADIENTS);
if(MYID==0) printf("\n This is an experimental feature.");
for (i=1;i<=NX;i++){
for (j=1;j<=NY;j++){
if(!ACOUSTIC) grad_vs[j][i]=grad_vs[j][i]*LBFGS_SCALE_GRADIENTS;
if(NPAR_LBFGS>1) grad_rho[j][i]=grad_rho[j][i]*LBFGS_SCALE_GRADIENTS;
if(NPAR_LBFGS>2) grad_vp[j][i]=grad_vp[j][i]*LBFGS_SCALE_GRADIENTS;
}
}
}
/*-------------------------------------------------*/
/* Init L-BFGS at iter==LBFGS_iter_start */
/*-------------------------------------------------*/
......
src/calc_energy.c
View file @ 0904becf
......@@ -30,8 +30,8 @@ double calc_energy(float **sectiondata, int ntr, int ns, float energy, int ntr_g
extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY;
int i, j;
float energy_dummy, intseis;
int umax=0, h;
float energy_dummy;
int h;
float intseis_sd;
float **intseis_sectiondata=NULL;
float *picked_times=NULL;
......@@ -59,11 +59,19 @@ double calc_energy(float **sectiondata, int ntr, int ns, float energy, int ntr_g
kill_tmp = imatrix(1,ntr_glob,1,nsrc_glob);
kill_vector = ivector(1,ntr);
if(TRKILL_OFFSET) {
/* Generate TraceKill file on the fly */
/*------------------*/
/* clear kill table */
/*------------------*/
for(i=1;i<=nsrc_glob;i++){
for (j=1; j<=ntr_glob; j++) {
kill_tmp[j][i]=0;
}
}
/* Use ONLY offset based TraceKill */
if(TRKILL_OFFSET==1) {
/* Generate TraceKill file on the fly based on the offset from the source */
create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,ishot,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
} else {
/* READ TraceKill file from disk */
......@@ -72,6 +80,7 @@ double calc_energy(float **sectiondata, int ntr, int ns, float energy, int ntr_g
sprintf(trace_kill_file,"%s_%i.dat",TRKILL_FILE,WORKFLOW_STAGE);
ftracekill=fopen(trace_kill_file,"r");
if (ftracekill==NULL){
/* If Workflow TraceKill file not found use File without workflow extensions */
sprintf(trace_kill_file,"%s.dat",TRKILL_FILE);
ftracekill=fopen(trace_kill_file,"r");
if (ftracekill==NULL){
......@@ -89,11 +98,18 @@ double calc_energy(float **sectiondata, int ntr, int ns, float energy, int ntr_g
for(i=1;i<=ntr_glob;i++){
for(j=1;j<=nsrc_glob;j++){
fscanf(ftracekill,"%d",&kill_tmp[i][j]);
if(feof(ftracekill)){
declare_error(" Error while reading TraceKill file. Check dimensions!");
}
}
}
fclose(ftracekill);
/* Use Tracekill FILE and add the offset based TraceKill */
if(TRKILL_OFFSET==2) {
create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,ishot,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
}
}
/* Generate local kill vector */
......
src/calc_hilbert.c
View file @ 0904becf
......@@ -25,8 +25,8 @@
void calc_hilbert(float ** datatrace, float ** envelope, int ns, int ntr){
/* declaration of variables */
int i,j, nfreq, npad, k;
float xr, yr, x, y, dump, a, *h;
int i,j, npad, k;
float *h;
double npadd;
/* declaration of variables for FFTW3*/
......
src/calc_misfit.c
View file @ 0904becf
......@@ -53,7 +53,6 @@ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LN
extern int TRKILL_OFFSET;
extern float TRKILL_OFFSET_LOWER;
extern float TRKILL_OFFSET_UPPER;
/*-----------------*/
/* Start Tracekill */
/*-----------------*/
......@@ -61,19 +60,28 @@ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LN
kill_tmp = imatrix(1,ntr_glob,1,nsrc_glob);
kill_vector = ivector(1,ntr);
if(TRKILL_OFFSET) {
/* Generate TraceKill file on the fly */
/*------------------*/
/* clear kill table */
/*------------------*/
for(i=1;i<=nsrc_glob;i++){
for (j=1; j<=ntr_glob; j++) {
kill_tmp[j][i]=0;
}
}
/* Use ONLY offset based TraceKill */
if(TRKILL_OFFSET==1) {
/* Generate TraceKill file on the fly based on the offset from the source */
create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,ishot,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
} else {
/* READ TraceKill file from disk */
if(USE_WORKFLOW){
sprintf(trace_kill_file,"%s_%i.dat",TRKILL_FILE,WORKFLOW_STAGE);
ftracekill=fopen(trace_kill_file,"r");
if (ftracekill==NULL){
/* If Workflow TraceKill file not found use File without workflow extensions */
sprintf(trace_kill_file,"%s.dat",TRKILL_FILE);
ftracekill=fopen(trace_kill_file,"r");
if (ftracekill==NULL){
......@@ -91,11 +99,18 @@ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LN
for(i=1;i<=ntr_glob;i++){
for(j=1;j<=nsrc_glob;j++){
fscanf(ftracekill,"%d",&kill_tmp[i][j]);
if(feof(ftracekill)){
declare_error(" Error while reading TraceKill file. Check dimensions!");
}
}
}
fclose(ftracekill);
/* Use Tracekill FILE and add the offset based TraceKill */
if(TRKILL_OFFSET==2) {
create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,-100,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
}
}
/* Generate local kill vector */
......
src/calc_res.c
View file @ 0904becf
......@@ -29,8 +29,8 @@ double calc_res(float **sectiondata, float **section, float **sectiondiff, float
extern int TRKILL, NORMALIZE, FC, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY, USE_WORKFLOW, WORKFLOW_STAGE;
float RMS, signL1, intseis;
int Lcount,i,j,invtime,k,h, umax=0;
float RMS, signL1;
int i,j,invtime,h, umax=0;
float l2;
float abs_section, abs_sectiondata, sectiondata_mult_section;
float intseis_s, intseis_sd;
......@@ -69,11 +69,19 @@ double calc_res(float **sectiondata, float **section, float **sectiondiff, float
kill_tmp = imatrix(1,ntr_glob,1,nsrc_glob);
kill_vector = ivector(1,ntr);
if(TRKILL_OFFSET) {
/* Generate TraceKill file on the fly */
/*------------------*/
/* clear kill table */
/*------------------*/
for(i=1;i<=nsrc_glob;i++){
for (j=1; j<=ntr_glob; j++) {
kill_tmp[j][i]=0;
}
}
/* Use ONLY offset based TraceKill */
if(TRKILL_OFFSET==1) {
/* Generate TraceKill file on the fly based on the offset from the source */
create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,ishot,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
} else {
/* READ TraceKill file from disk */
......@@ -82,6 +90,7 @@ double calc_res(float **sectiondata, float **section, float **sectiondiff, float
sprintf(trace_kill_file,"%s_%i.dat",TRKILL_FILE,WORKFLOW_STAGE);
ftracekill=fopen(trace_kill_file,"r");
if (ftracekill==NULL){
/* If Workflow TraceKill file not found use File without workflow extensions */
sprintf(trace_kill_file,"%s.dat",TRKILL_FILE);
ftracekill=fopen(trace_kill_file,"r");
if (ftracekill==NULL){
......@@ -99,11 +108,18 @@ double calc_res(float **sectiondata, float **section, float **sectiondiff, float
for(i=1;i<=ntr_glob;i++){
for(j=1;j<=nsrc_glob;j++){
fscanf(ftracekill,"%d",&kill_tmp[i][j]);
if(feof(ftracekill)){
declare_error(" Error while reading TraceKill file. Check dimensions!");
}
}
}
fclose(ftracekill);
/* Use Tracekill FILE and add the offset based TraceKill */
if(TRKILL_OFFSET==2) {
create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,-100,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
}
}
/* Generate local kill vector */
......@@ -117,9 +133,6 @@ double calc_res(float **sectiondata, float **section, float **sectiondiff, float
/* End Tracekill */
/*---------------*/
RMS=0.0;
Lcount=1;
/* Integration of measured and synthetic data */
for(i=1;i<=ntr;i++){
intseis_s=0;
......
src/count_killed_traces.c
View file @ 0904becf
......@@ -39,20 +39,28 @@ void count_killed_traces(int ntr, int swstestshot, int ntr_glob, int **recpos_lo
FILE *ftracekill;
kill_tmp = imatrix(1,ntr_glob,1,nsrc_glob);
kill_vector = ivector(1,ntr);
/*------------------*/
/* clear kill table */
/*------------------*/
for(i=1;i<=nsrc_glob;i++){
for (j=1; j<=ntr_glob; j++) {
kill_tmp[j][i]=0;
}
}
if(TRKILL_OFFSET) {
/* Use ONLY offset based TraceKill */
if(TRKILL_OFFSET==1) {
/* Generate TraceKill file on the fly based on the offset from the source */
create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,ishot,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
if(MYID==0) {
printf("\n\n ----- Offset based Tracekill ------");
printf("\n Kill offsets between %.1f m and %.1f m",TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
}
/* Generate TraceKill file on the fly */
create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,ishot,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
} else {
/* READ TraceKill file from disk */
......@@ -61,6 +69,7 @@ void count_killed_traces(int ntr, int swstestshot, int ntr_glob, int **recpos_lo
sprintf(trace_kill_file,"%s_%i.dat",TRKILL_FILE,WORKFLOW_STAGE);
ftracekill=fopen(trace_kill_file,"r");
if (ftracekill==NULL){
/* If Workflow TraceKill file not found use File without workflow extensions */
sprintf(trace_kill_file,"%s.dat",TRKILL_FILE);
ftracekill=fopen(trace_kill_file,"r");
if (ftracekill==NULL){
......@@ -78,14 +87,26 @@ void count_killed_traces(int ntr, int swstestshot, int ntr_glob, int **recpos_lo
for(i=1;i<=ntr_glob;i++){
for(j=1;j<=nsrc_glob;j++){
fscanf(ftracekill,"%d",&kill_tmp[i][j]);
if(feof(ftracekill)){
declare_error(" Error while reading TraceKill file. Check dimensions!");
}
}
}
fclose(ftracekill);
/* Use Tracekill FILE and add the offset based TraceKill */
if(TRKILL_OFFSET==2) {
create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,-100,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
if(MYID==0) {
printf("\n\n ----- Offset based Tracekill ------");
printf("\n Kill offsets between %.1f m and %.1f m",TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
printf("\n In addition the TraceKill File is used");
}
}
}
h=1;
for(i=1;i<=ntr;i++){
......
src/create_trkill_table.c
View file @ 0904becf
......@@ -27,15 +27,7 @@ void create_trkill_table(int ** killtable, int ntr_glob, int **recpos, int nsrc_
/* extern variables */
extern float DH;
extern int MYID;
/*------------------*/
/* clear kill table */
/*------------------*/
for(s=1;s>=nsrc_glob;s++){
for (r=1; r<=ntr_glob; r++) {
killtable[r][s]=0;
}
}
extern char TRKILL_FILE[STRING_SIZE];
/*---------------------*/
/* Generate Killtable */
......@@ -57,11 +49,11 @@ void create_trkill_table(int ** killtable, int ntr_glob, int **recpos, int nsrc_
/* Debug */
/*---------------------*/
/* Debug is declared with ishot=-100 */
if(ishot==-100){
if((ishot==-100) && (MYID==0)){
FILE *FP_TK_OUT;
char filename[225];
sprintf(filename,"tracekill_%.0f_%.0f.txt",kill_offset_lower,kill_offset_upper);
sprintf(filename,"%stracekill.out.%.0f_%.0f.txt",TRKILL_FILE,kill_offset_lower,kill_offset_upper);
FP_TK_OUT=fopen(filename,"w");
for (r=1; r<=ntr_glob; r++) {
......
src/exchange_par.c
View file @ 0904becf
......@@ -107,6 +107,7 @@ void exchange_par(void){
extern int EPRECOND_PER_SHOT;
extern int EPRECOND_PER_SHOT_SH;
extern float LBFGS_SCALE_GRADIENTS;
extern int LBFGS_SURFACE;
extern int LBFGS_STEP_LENGTH;
......@@ -217,6 +218,8 @@ void exchange_par(void){
fdum[67]=TRKILL_STF_OFFSET_UPPER;
fdum[68]=TRKILL_OFFSET_LOWER;
fdum[69]=TRKILL_OFFSET_UPPER;
fdum[70]=LBFGS_SCALE_GRADIENTS;
/***********/
/* Integer */
......@@ -499,6 +502,8 @@ void exchange_par(void){
TRKILL_OFFSET_LOWER=fdum[68];
TRKILL_OFFSET_UPPER=fdum[69];
LBFGS_SCALE_GRADIENTS=fdum[70];
/***********/
/* Integer */
/***********/
......
src/globvar.h
View file @ 0904becf
......@@ -93,6 +93,7 @@ int N_LBFGS;
int LBFGS_SURFACE;
int LBFGS_STEP_LENGTH;
float LBFGS_SCALE_GRADIENTS;
int WOLFE_CONDITION=0;
int WOLFE_NUM_TEST=5;
int WOLFE_TRY_OLD_STEPLENGTH=1;
......
src/read_par_json.c
View file @ 0904becf
......@@ -121,6 +121,7 @@ void read_par_json(FILE *fp, char *fileinp){
extern int LBFGS_STEP_LENGTH;
extern int N_LBFGS;
extern float LBFGS_SCALE_GRADIENTS;
extern int WOLFE_CONDITION;
extern int WOLFE_NUM_TEST;
extern int WOLFE_TRY_OLD_STEPLENGTH;
......@@ -712,6 +713,10 @@ void read_par_json(FILE *fp, char *fileinp){
fprintf(fp,"Variable N_LBFGS is set to default value %d.\n",N_LBFGS);
}
if (get_float_from_objectlist("LBFGS_SCALE_GRADIENTS",number_readobjects,&LBFGS_SCALE_GRADIENTS,varname_list, value_list)){
LBFGS_SCALE_GRADIENTS=1;
}
if (get_int_from_objectlist("WOLFE_CONDITION",number_readobjects,&WOLFE_CONDITION,varname_list, value_list)){
WOLFE_CONDITION=1;
fprintf(fp,"Variable WOLFE_CONDITION is set to default value %d.\n",WOLFE_CONDITION);
......@@ -765,6 +770,10 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_float_from_objectlist("TRKILL_STF_OFFSET_UPPER",number_readobjects,&TRKILL_STF_OFFSET_UPPER,varname_list, value_list)){
declare_error("Variable TRKILL_STF_OFFSET_UPPER could not be retrieved from the json input file!");
}
if(TRKILL_STF_OFFSET==2){
if (get_string_from_objectlist("TRKILL_FILE_STF",number_readobjects,TRKILL_FILE_STF,varname_list, value_list))
declare_error("Variable TRKILL_FILE_STF could not be retrieved from the json input file!");
}
}
}
}
......@@ -848,7 +857,6 @@ void read_par_json(FILE *fp, char *fileinp){
fprintf(fp,"Variable TRKILL is set to default value %d.\n",TRKILL);}
else {
if (TRKILL==1) {
if (get_int_from_objectlist("TRKILL_OFFSET",number_readobjects,&TRKILL_OFFSET,varname_list, value_list)){
TRKILL_OFFSET=0;
if (get_string_from_objectlist("TRKILL_FILE",number_readobjects,TRKILL_FILE,varname_list, value_list))
......@@ -860,6 +868,10 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_float_from_objectlist("TRKILL_OFFSET_UPPER",number_readobjects,&TRKILL_OFFSET_UPPER,varname_list, value_list)){
declare_error("Variable TRKILL_OFFSET_UPPER could not be retrieved from the json input file!");
}
if(TRKILL_OFFSET==2){
if (get_string_from_objectlist("TRKILL_FILE",number_readobjects,TRKILL_FILE,varname_list, value_list))
declare_error("Variable TRKILL_FILE could not be retrieved from the json input file!");
}
}
}
}
......
src/stf.c
View file @ 0904becf
......@@ -70,7 +70,17 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy
kill_tmp = imatrix(1,ntr_glob,1,nshots);
kill_vector = ivector(1,ntr_glob);
if(TRKILL_STF_OFFSET) {
/*------------------*/
/* clear kill table */
/*------------------*/
for(i=1;i<=nsrc_glob;i++){
for (j=1; j<=ntr_glob; j++) {
kill_tmp[j][i]=0;
}
}
if(TRKILL_STF_OFFSET==1) {
if(MYID==0) {
printf("Automatic offset based TraceKill for STF\n");
......@@ -105,6 +115,11 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy
}
fclose(ftracekill);
if(TRKILL_STF_OFFSET==2) {
/* Generate TraceKill file on the fly */
create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,-100,TRKILL_STF_OFFSET_LOWER,TRKILL_STF_OFFSET_UPPER);
}
}
h=1;
for(i=1;i<=ntr_glob;i++){
......
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