Offset based tracekill STF + Makefile improvement

Implemented an offset based tracekill for STF.
Makefile in par folder has now two new targets:
1. make install, do an install from precompiled files
2. make reinstall, do a clean install

genmod/toy_example_elastic_manuell.c

+/*-----------------------------------------------------------------------------------------
+ * Copyright (C) 2016
+ *
+ * This file is part of IFOS2D.
+ *
+ * IFOS2D is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, version 2.0 of the License only.
+ *
+ * IFOS2D is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with IFOS2D. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
+ -----------------------------------------------------------------------------------------*/
+
+
+#include "fd.h"
+
+void model_elastic(float  **  rho, float **  pi, float **  u){
+    
+    /*--------------------------------------------------------------------------*/
+    /* extern variables */
+    
+    extern int NX, NY, NXG, NYG,  POS[3], MYID, L;
+    extern float DH, DT, *FL, TAU;
+    extern char  MFILE[STRING_SIZE];
+    extern char INV_MODELFILE[STRING_SIZE];
+    
+    /* local variables */
+    float muv, piv, vp, vs, rhov, ts, tp, *pts;
+    int i, j, ii, jj, l;
+    char modfile[STRING_SIZE];
+    
+    FILE *flfile;
+    int nodes;
+    char cline[256];
+    
+    float *fldepth, *flrho, *flvp, *flvs;
+    
+    
+    
+    
+    /**************************************************/
+    /* creation of shear wave velocity and TAU models */
+    /**************************************************/
+    
+    /*read FL nodes from File*/
+    nodes=3;
+    fldepth=vector(1,nodes);
+    flvs=vector(1,nodes);
+    flrho=vector(1,nodes);
+    flvp=vector(1,nodes);
+    flfile=fopen("model_true/flnodes.toy_example.start","r");
+    if (flfile==NULL) declare_error(" FL-file could not be opened !");
+    
+    
+    
+    for (l=1;l<=nodes;l++){
+        fgets(cline,255,flfile);
+        if (cline[0]!='#'){
+            sscanf(cline,"%f%f%f%f",&fldepth[l], &flrho[l], &flvp[l], &flvs[l]);
+        }
+        else l=l-1;
+        
+    }
+    
+    if(MYID==0){
+        printf(" ------------------------------------------------------------------ \n\n");
+        printf(" Information of FL nodes: \n\n");
+        printf(" \t depth \t vs \n\n");
+        
+        for (l=1;l<=nodes;l++){
+            printf(" \t %f \t %f\n\n",fldepth[l],flvs[l]);
+        }
+        printf(" ------------------------------------------------------------------ \n\n");
+    }
+    
+    /* loop over global grid */
+    for (i=1;i<=NXG;i++){
+        for (l=1;l<nodes;l++){
+            if (fldepth[l]==fldepth[l+1]){
+                if ((i==1) && (MYID==0)){
+                    printf("depth: %f m: double node\n",fldepth[l]);}}
+            else{
+                for (j=(int)(fldepth[l]/DH)+1;j<=(int)(fldepth[l+1]/DH);j++){
+                    
+                    vs=0.0;
+                    
+                    vs=(DH*(j-1)-fldepth[l])*(flvs[l+1]-flvs[l])/(fldepth[l+1]-fldepth[l])+flvs[l];
+                    vs=vs*1000.0;
+                    
+                    muv=vs;
+                    ts=TAU;
+                    tp=TAU;
+                    
+                    /* only the PE which belongs to the current global gridpoint
+                     is saving model parameters in his local arrays */
+                    if ((POS[1]==((i-1)/NX)) &&
+                        (POS[2]==((j-1)/NY))){
+                        ii=i-POS[1]*NX;
+                        jj=j-POS[2]*NY;
+                        
+                        u[jj][ii]=muv;
+                    }
+                }
+            }
+        }
+        
+        for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
+            
+            vs=0.0;
+            vs=flvs[nodes]*1000.0;
+            
+            muv=vs;
+            
+            /* only the PE which belongs to the current global gridpoint
+             is saving model parameters in his local arrays */
+            if ((POS[1]==((i-1)/NX)) &&
+                (POS[2]==((j-1)/NY))){
+                ii=i-POS[1]*NX;
+                jj=j-POS[2]*NY;
+                
+                u[jj][ii]=muv;
+            }
+        }
+    }
+    
+    free_vector(fldepth,1,nodes);
+    free_vector(flrho,1,nodes);
+    free_vector(flvp,1,nodes);
+    free_vector(flvs,1,nodes);
+    
+    
+    
+    /**************************************************/
+    /* creation of P wave velocity and density models */
+    /**************************************************/
+    
+    /*read FL nodes from File*/
+    nodes=7;
+    fldepth=vector(1,nodes);
+    flvs=vector(1,nodes);
+    flrho=vector(1,nodes);
+    flvp=vector(1,nodes);
+    flfile=fopen("model_true/flnodes.toy_example","r");
+    flfile=fopen("model_true/flnodes.toy_example.start","r");
+    if (flfile==NULL) declare_error(" FL-file could not be opened !");
+    
+    for (l=1;l<=nodes;l++){
+        fgets(cline,255,flfile);
+        if (cline[0]!='#'){
+            sscanf(cline,"%f%f%f%f",&fldepth[l], &flrho[l], &flvp[l], &flvs[l]);
+        }
+        else l=l-1;
+    }
+    
+    
+    if(MYID==0){
+        printf(" ------------------------------------------------------------------ \n\n");
+        printf(" Information of FL nodes: \n\n");
+        printf(" \t depth \t vp \t rho \n\n");
+        
+        for (l=1;l<=nodes;l++){
+            printf(" \t %f \t %f \t %f \n\n",fldepth[l],flvp[l],flrho[l]);
+        }
+        printf(" ------------------------------------------------------------------ \n\n");
+    }
+    
+    /* loop over global grid */
+    for (i=1;i<=NXG;i++){
+        for (l=1;l<nodes;l++){
+            if (fldepth[l]==fldepth[l+1]){
+                if ((i==1) && (MYID==0)){
+                    printf("depth: %f m: double node\n",fldepth[l]);}}
+            else{
+                for (j=(int)(fldepth[l]/DH)+1;j<=(int)(fldepth[l+1]/DH);j++){
+                    
+                    vp=0.0; rhov=0.0;
+                    
+                    vp=(DH*(j-1)-fldepth[l])*(flvp[l+1]-flvp[l])/(fldepth[l+1]-fldepth[l])+flvp[l];
+                    vp=vp*1000.0;
+                    rhov=(DH*(j-1)-fldepth[l])*(flrho[l+1]-flrho[l])/(fldepth[l+1]-fldepth[l])+flrho[l];
+                    rhov=rhov*1000.0;
+                    
+                    piv=vp;
+                    
+                    /* only the PE which belongs to the current global gridpoint
+                     is saving model parameters in his local arrays */
+                    if ((POS[1]==((i-1)/NX)) &&
+                        (POS[2]==((j-1)/NY))){
+                        ii=i-POS[1]*NX;
+                        jj=j-POS[2]*NY;
+                        
+                        rho[jj][ii]=rhov;
+                        pi[jj][ii]=piv;
+                    }
+                }
+            }
+        }
+        
+        for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
+            
+            vp=0.0; rhov=0.0;
+            vp=flvp[nodes]*1000.0; rhov=flrho[nodes]*1000.0;
+            
+            piv=vp;
+            
+            /* only the PE which belongs to the current global gridpoint 
+             is saving model parameters in his local arrays */
+            if ((POS[1]==((i-1)/NX)) && 
+                (POS[2]==((j-1)/NY))){
+                ii=i-POS[1]*NX;
+                jj=j-POS[2]*NY;
+                
+                rho[jj][ii]=rhov;
+                pi[jj][ii]=piv;
+            }
+        }
+    }
+    
+//    free_vector(fldepth,1,nodes);
+//    free_vector(flrho,1,nodes);
+//    free_vector(flvp,1,nodes);
+//    free_vector(flvs,1,nodes);
+//    
+    
+    
+    
+    sprintf(modfile,"%s_rho_it0.bin",INV_MODELFILE);
+    writemod(modfile,rho,3);
+    MPI_Barrier(MPI_COMM_WORLD);
+    if (MYID==0) mergemod(modfile,3);
+    MPI_Barrier(MPI_COMM_WORLD); 
+    sprintf(modfile,"%s_rho_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
+    remove(modfile);
+    
+    sprintf(modfile,"%s_vs_it0.bin",INV_MODELFILE);
+    writemod(modfile,u,3);
+    MPI_Barrier(MPI_COMM_WORLD);
+    if (MYID==0) mergemod(modfile,3);
+    MPI_Barrier(MPI_COMM_WORLD); 
+    sprintf(modfile,"%s_vs_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
+    remove(modfile);
+    
+    sprintf(modfile,"%s_vp_it0.bin",INV_MODELFILE);
+    writemod(modfile,pi,3);
+    MPI_Barrier(MPI_COMM_WORLD);
+    if (MYID==0) mergemod(modfile,3);
+    MPI_Barrier(MPI_COMM_WORLD); 
+    sprintf(modfile,"%s_vp_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
+    remove(modfile);
+    
+    
+    
+    
+//    free_vector(pts,1,L);
+    
+}
+
+
+

par/Makefile

@@ -26,7 +26,7 @@ else
 endif
 
 .PHONY: install
-install: clean IFOS2D
+install: clean-bin IFOS2D
 
 .PHONY: clean
 clean: 
@@ -36,3 +36,10 @@ clean:
 	$(MAKE) -C ../contrib/stfinv clean
 	$(MAKE) -C ../contrib/libcseife clean
 	$(MAKE) -C ../src/ clean
+
+.PHONY: clean-bin
+clean-bin:
+	rm -rf ../bin/*
+
+.PHONY: reinstall
+reinstall: clean IFOS2D
\ No newline at end of file

par/in_and_out/IFOS2D_INV_all_parameters.json

@@ -155,6 +155,8 @@
 			"N_STF_START" : "1",
 			"TAPER_STF" : "0",
 			"TRKILL_STF" : "0",
+			"TRKILL_STF_OFFSET" : "0",
+			"TRKILL_STF_OFFSET_LOWER , TRKILL_STF_OFFSET_UPPER" : "0 , 10",
 			"TRKILL_FILE_STF" : "./trace_kill/trace_kill.dat",
 
 "Frequency filtering during inversion" : "comment",
@@ -181,6 +183,8 @@
 
 "Trace killing" : "comment",
 			"TRKILL" : "0",
+			"TRKILL_OFFSET" : "0",
+			"TRKILL_OFFSET_LOWER , TRKILL_OFFSET_UPPER" : "0 , 10",
 			"TRKILL_FILE" : "./trace_kill/trace_kill",
 
 "Definition of a gradient taper" : "comment",

par/run_toy_example.sh

@@ -5,8 +5,7 @@
 ###############################################################
 
 # compiling all libraries and  IFOS
-make clean
-make IFOS2D MODEL=../genmod/toy_example_true.c
+make install MODEL=../genmod/toy_example_true.c
 
 # starting IFOS for forward modeling
 # (depending on the MPI package you maybe have to adjust the following programme call, 
@@ -19,8 +18,7 @@ mpirun -np 4 nice -19 ../bin/IFOS2D in_and_out/toy_example/toy_example_FW.json |
 ###############################################################
 
 # compiling IFOS
-make clean
-make IFOS2D MODEL=../genmod/toy_example_start.c
+make install MODEL=../genmod/toy_example_start.c
 
 # starting IFOS
 #lamboot

par/run_toy_example_SH.sh

@@ -5,8 +5,7 @@
 ###############################################################
 
 # compiling all libraries and  IFOS
-make clean
-make IFOS2D MODEL=../genmod/toy_example_true.c MODEL_EL=../genmod/toy_example_elastic_true.c
+make install MODEL=../genmod/toy_example_true.c MODEL_EL=../genmod/toy_example_elastic_true.c
 
 # starting IFOS for forward modeling
 # (depending on the MPI package you maybe have to adjust the following programme call, 
@@ -20,8 +19,7 @@ mpirun -np 4 ../bin/IFOS2D in_and_out/toy_example/toy_example_FW_SH.json | tee i
 ###############################################################
 
 # compiling IFOS
-make clean
-make IFOS2D MODEL=../genmod/toy_example_start.c MODEL_EL=../genmod/toy_example_elastic_start.c
+make install MODEL=../genmod/toy_example_start.c MODEL_EL=../genmod/toy_example_elastic_start.c
 
 # starting IFOS
 #lamboot

par/run_toy_example_ac.sh

@@ -12,8 +12,7 @@ rm su/toy_example/toy_example*
 rm su/measured_data/toy_example*
 
 # compiling all libraries and  IFOS
-make clean
-make IFOS2D MODEL_AC=../genmod/toy_example_ac_true.c
+make install MODEL_AC=../genmod/toy_example_ac_true.c
 
 # starting IFOS for forward modeling
 mpirun -np 4 nice -19 ../bin/IFOS2D in_and_out/toy_example/toy_example_ac_FW.json | tee in_and_out/toy_example/toy_example_ac_FW.out
@@ -23,8 +22,7 @@ mpirun -np 4 nice -19 ../bin/IFOS2D in_and_out/toy_example/toy_example_ac_FW.jso
 ###############################################################
 
 # compiling IFOS
-make clean
-make IFOS2D MODEL_AC=../genmod/toy_example_ac_start.c
+make install MODEL_AC=../genmod/toy_example_ac_start.c
 
 # starting IFOS
 mpirun -np 4 nice -19 ../bin/IFOS2D in_and_out/toy_example/toy_example_ac_INV.json | tee in_and_out/toy_example/toy_example_ac_INV.out

src/IFOS2D.c

@@ -32,7 +32,7 @@
 
 int main(int argc, char **argv){
     /* variables in main */
-    int ns, nseismograms=0, nt, nd, fdo3, j, i, iter, h, infoout, SHOTINC,  hin, hin1, s=0;
+    int ns, nseismograms=0, nt, nd, fdo3, j, i, iter, h, infoout, SHOTINC,  hin, hin1, do_stf=0;
     int NTDTINV, nxny, nxnyi, imat, imat1, imat2, IDXI, IDYI, hi, NTST, NTSTI;
     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, itestshot;
@@ -1172,7 +1172,7 @@ int main(int argc, char **argv){
                          * Therefore (gradient_optimization==1) is added.
                          */
                         
-                        if(((INV_STF==1)&&((iter==1)||(s==1))) && (gradient_optimization==1)){
+                        if(((INV_STF==1)&&( (iter==1) || (do_stf==1) )) && (gradient_optimization==1) ){
                             fprintf(FP,"\n==================================================================================\n");
                             fprintf(FP,"\n MYID=%d *****  Forward simulation for inversion of source time function ******** \n",MYID);
                             fprintf(FP,"\n MYID=%d * Starting simulation (forward model) for shot %d of %d. Iteration %d ** \n",MYID,ishot,nshots,iter);
@@ -1441,7 +1441,7 @@ int main(int argc, char **argv){
                             if((TIME_FILT==1) ||(TIME_FILT==2)){
                                 
                                 if (!FORWARD_ONLY){
-                                    if((INV_STF==1)&&((iter==1)||(s==1))){
+                                    if((INV_STF==1)&&((iter==1)||(do_stf==1))){
                                         
                                         if (nsrc_loc>0){
                                             
@@ -1477,21 +1477,26 @@ int main(int argc, char **argv){
                                             
                                             if(WAVETYPE==1 || WAVETYPE==3){
                                                 if ((ADJOINT_TYPE==1)|| (ADJOINT_TYPE==2)){
-                                                    stf(FP,fulldata_vy,sectionvy_obs,sectionvy_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0);
+                                                    stf(FP,fulldata_vy,sectionvy_obs,sectionvy_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0,nsrc_glob);
                                                 }
                                                 if (ADJOINT_TYPE==4){
-                                                    stf(FP,fulldata_p,sectionp_obs,sectionp_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0);
+                                                    stf(FP,fulldata_p,sectionp_obs,sectionp_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0,nsrc_glob);
                                                 }
                                             }
                                             
                                             if(WAVETYPE==2 || WAVETYPE==3){
-                                                stf(FP,fulldata_vz,sectionvz_obs,sectionvz_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,1);
+                                                stf(FP,fulldata_vz,sectionvz_obs,sectionvz_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,1,nsrc_glob);
                                             }
+                                            
+                                            /* do_stf controles if a STF is done or not */
+                                            do_stf=0;
+                                            /* As STF was done no, it is deactivated (do_stf=0) until it will be activated*/
                                         }
                                     }
                                 }
                             
                             } else {
+                                
                                 if (FORWARD_ONLY==0){
                                     if((INV_STF==1)&&(iter==N_STF_START)){
                                         
@@ -1512,15 +1517,15 @@ int main(int argc, char **argv){
                                                 }
                                                 
                                                 if ((ADJOINT_TYPE==1)|| (ADJOINT_TYPE==2)){
-                                                    stf(FP,fulldata_vy,sectionvy_obs,sectionvy_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0);
+                                                    stf(FP,fulldata_vy,sectionvy_obs,sectionvy_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0,nsrc_glob);
                                                 }
                                                 if (ADJOINT_TYPE==4){
-                                                    stf(FP,fulldata_p,sectionp_obs,sectionp_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0);
+                                                    stf(FP,fulldata_p,sectionp_obs,sectionp_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,0,nsrc_glob);
                                                 }
                                             }
                                             if(WAVETYPE==2 || WAVETYPE==3){
                                                 inseis(fprec,ishot,sectionvz_obs,ntr_glob,ns,10,iter);
-                                                stf(FP,fulldata_vz,sectionvz_obs,sectionvz_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,1);
+                                                stf(FP,fulldata_vz,sectionvz_obs,sectionvz_conv,source_time_function,recpos,recpos_loc,ntr_glob,ntr,srcpos,ishot,ns,iter,nsrc_glob,FC,1,nsrc_glob);
                                             }
                                         }
                                     }
@@ -3760,7 +3765,7 @@ int main(int argc, char **argv){
             
             /* saving history of final L2*/
             L2_hist[iter]=L2t[4];
-            s=0;
+            do_stf=0;
             
             
             /* -----------------------------------------------------------------------*/
@@ -3869,7 +3874,7 @@ int main(int argc, char **argv){
                 fprintf(FP,"\n Switching to next line in workflow");
                 
                 WORKFLOW_STAGE++;
-                s=1;
+                do_stf=1;
                 min_iter_help=0;
                 min_iter_help=iter+MIN_ITER;
                 
@@ -3899,7 +3904,7 @@ int main(int argc, char **argv){
                 }
                 
                 FC=FC+FC_INCR;
-                s=1;
+                do_stf=1;
                 min_iter_help=0;
                 min_iter_help=iter+MIN_ITER;
                 
@@ -3938,7 +3943,7 @@ int main(int argc, char **argv){
                 
                 FREQ_NR=FREQ_NR+1;
                 FC=FC_EXT[FREQ_NR];
-                s=1;
+                do_stf=1;
                 min_iter_help=0;
                 min_iter_help=iter+MIN_ITER;
                 if(MYID==0) printf("\n Changing to corner frequency of %4.2f Hz \n",FC);

src/calc_energy.c

@@ -60,11 +60,7 @@ double calc_energy(float **sectiondata, int ntr, int ns, float energy, int ntr_g
         kill_vector = ivector(1,ntr);
         
         if(TRKILL_OFFSET) {
-            
-            if(MYID==0) {
-                printf("Automatic offset based TraceKill: Experimental feature");
-            }
-            
+
             /* Generate TraceKill file on the fly */
             create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,ishot,TRKILL_OFFSET_LOWER,TRKILL_OFFSET_UPPER);
             

src/calc_misfit.c

@@ -63,10 +63,6 @@ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LN
         
         if(TRKILL_OFFSET) {
             
-            if(MYID==0) {
-                printf("Automatic offset based TraceKill: Experimental feature");
-            }
-            
             /* Generate TraceKill file on the fly */