JOINT: EPRECOND & joint function

Separate waterlevel for SH-waves && new joint function which takes the mean without normalization.

src/denise.c

@@ -2337,7 +2337,6 @@ int main(int argc, char **argv){
 
                                         hin++;
                                     }
-                                    /*fclose(FP2);*/
                                     
                                     if((EPRECOND==1)&&(EPRECOND_ITER==iter||(EPRECOND_ITER==0))){
                                         if (WAVETYPE==1 || WAVETYPE==3) eprecond(Wr,pvx,pvy);
@@ -2551,7 +2550,6 @@ int main(int argc, char **argv){
                                         
                                         if(INVMAT1==1){
                                             /* calculate density gradient */
-                                            
                                             waveconv_rho_shot_z[j][i] = ( (pu[j][i] * pu[j][i] * waveconv_mu_z[j][i]) + waveconv_rho_s_z[j][i]);
                                             
                                         }
@@ -2570,15 +2568,18 @@ int main(int argc, char **argv){
                             /* calculate and apply energy preconditioning   */
                             /* -------------------------------------------- */
                             if((EPRECOND==1)||(EPRECOND==3)){
+                                
                                 /* calculate energy weights */
                                 if(EPRECOND_ITER==iter||(EPRECOND_ITER==0)) {
                                     fprintf(FP,"\n Calculating approx. Hessian for shot %i. EPRECOND=%i, EPSILON_WE=%f",ishot,EPRECOND,EPSILON_WE);
+                                    
                                     if(WAVETYPE==1 || WAVETYPE==3) {
-                                        eprecond1(We,Ws,Wr);
+                                        eprecond1(We,Ws,Wr,EPSILON_WE);
                                         if(EPRECOND_PER_SHOT) We_max=global_maximum(We);
                                     }
+                                    
                                     if(WAVETYPE==2 || WAVETYPE==3) {
-                                        eprecond1(We_SH,Ws_SH,Wr_SH);
+                                        eprecond1(We_SH,Ws_SH,Wr_SH,EPSILON_WE_SH);
                                         if(EPRECOND_PER_SHOT) We_max_SH=global_maximum(We_SH);
                                     }
                                     

src/eprecond1.c

 /* This function is copied from DENISE Black Edition from D. Koehn  */
 
 #include "fd.h"
-void eprecond1(float ** We, float ** Ws, float ** Wr){
+void eprecond1(float ** We, float ** Ws, float ** Wr, float epsilon){
     
     extern int NX, NY, IDX, IDY, DTINV, EPRECOND, VERBOSE;
     extern int POS[3], NXG;
     extern float DH;
     int i, j, k, l, ii, jj;
     float maxWetmp, maxWe, x, y, xmin, xmax;
-    extern float EPSILON_WE;
     xmin = DH;
     xmax = NXG*DH;
     
@@ -56,11 +55,12 @@ void eprecond1(float ** We, float ** Ws, float ** Wr){
     
     /* estimate maximum of We */
     MPI_Allreduce(&maxWetmp,&maxWe,1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);
+    
     /* regularize energy weighting to avoid divison by zero */
     for (j=1;j<=NY;j=j+IDY){
         for (i=1;i<=NX;i=i+IDX){
             
-            We[j][i] = We[j][i] + (EPSILON_WE*maxWe);
+            We[j][i] = We[j][i] + (epsilon*maxWe);
             
         }
     }

src/exchange_par.c

@@ -100,7 +100,7 @@ void exchange_par(void){
     
     extern int EPRECOND;
     extern int EPRECOND_ITER;
-    extern float EPSILON_WE;
+    extern float EPSILON_WE,EPSILON_WE_SH;
     extern int EPRECOND_PER_SHOT;
     
     extern int LBFGS_SURFACE;
@@ -114,16 +114,18 @@ void exchange_par(void){
     extern float WOLFE_C1_SL;
     extern float WOLFE_C2_SL;
     
-    /* definition of local variables */
     
+    /* definition of local variables */
     /* NPAR is set in fd.h and must be set to the max. number of elements in idum/fdum +1 (because vector index starts at 0) */
-    
     int idum[NPAR];
     float fdum[NPAR];
     
     
     if (MYID == 0){
         
+        /***********/
+        /*  Float  */
+        /***********/
         fdum[1]  = DH;
         fdum[2]  = TIME;
         fdum[3]  = DT;
@@ -202,10 +204,14 @@ void exchange_par(void){
         fdum[61] = JOINT_INVERSION_PSV_SH_ALPHA_RHO;
         
         fdum[62]=EPSILON_WE;
+        fdum[63]=EPSILON_WE_SH;
         
-        fdum[63]=WOLFE_C1_SL;
-        fdum[64]=WOLFE_C2_SL;
+        fdum[64]=WOLFE_C1_SL;
+        fdum[65]=WOLFE_C2_SL;
         
+        /***********/
+        /* Integer */
+        /***********/
         idum[1]  = NPROCX;
         idum[2]  = NPROCY;
         
@@ -396,6 +402,9 @@ void exchange_par(void){
     
     MPI_Barrier(MPI_COMM_WORLD);
     
+    /***********/
+    /*  Float  */
+    /***********/
     DH=fdum[1];
     TIME=fdum[2];
     DT=fdum[3];
@@ -473,9 +482,14 @@ void exchange_par(void){
     JOINT_INVERSION_PSV_SH_ALPHA_RHO = fdum[61];
     
     EPSILON_WE=fdum[62];
+    EPSILON_WE_SH=fdum[63];
+    
+    WOLFE_C1_SL=fdum[64];
+    WOLFE_C2_SL=fdum[65];
     
-    WOLFE_C1_SL=fdum[63];
-    WOLFE_C2_SL=fdum[64];
+    /***********/
+    /* Integer */
+    /***********/
     
     NPROCX = idum[1];
     NPROCY = idum[2];

src/fd.h

@@ -527,7 +527,7 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST
 
 void eprecond(float ** W, float ** vx, float ** vy);
 void eprecond_SH(float ** W, float ** vz);
-void eprecond1(float ** We, float ** Ws, float ** Wr);
+void eprecond1(float ** We, float ** Ws, float ** Wr, float epsilon);
 
 /* Matrix Operations */
 float average_matrix(float ** matrix);

src/globvar.h

@@ -85,7 +85,7 @@ float EPS_SCALE, SCALEFAC;
 float PRO;
 
 int EPRECOND;
-float EPSILON_WE;
+float EPSILON_WE,EPSILON_WE_SH;
 int EPRECOND_ITER;
 int EPRECOND_PER_SHOT;
 

src/joint_inversion.c

@@ -25,9 +25,6 @@
  *  ----------------------------------------------------------------------*/
 #include "fd.h"
 
-/* Declare local functions */
-
-
 float ** joint_inversion_grad ( float ** gradiant_1,float ** gradiant_2, float alpha, int joint_type) {
     /*
      * Input parameters:
@@ -57,11 +54,31 @@ float ** joint_inversion_grad ( float ** gradiant_1,float ** gradiant_2, float a
             for (j=1;j<=NY;j++){
                 for (i=1;i<=NX;i++){
                     gradiant_1[j][i]=((1-alpha)*gradiant_1[j][i]/max1+alpha*gradiant_2[j][i]/max2);
+                    
+                    if (isnan(gradiant_1[j][i])) {
+                        gradiant_1[j][i]=0.0;
+                    }
                 }
             }
             joint_gradiant=gradiant_1;
             break;
+            
+        case 2:
+             /* joint_type=2: Arithmetic mean without normalization */
+            for (j=1;j<=NY;j++){
+                for (i=1;i<=NX;i++){
+                    gradiant_1[j][i]=((1-alpha)*gradiant_1[j][i]+alpha*gradiant_2[j][i]);
+                    
+                    if (isnan(gradiant_1[j][i])) {
+                        gradiant_1[j][i]=0.0;
+                    }
+                }
+            }
+            joint_gradiant=gradiant_1;
+            break;
+            
     }
+
     return joint_gradiant;
 };
 

src/read_par_json.c

@@ -107,7 +107,7 @@ void read_par_json(FILE *fp, char *fileinp){
     
     extern int EPRECOND;
     extern int EPRECOND_ITER;
-    extern float EPSILON_WE;
+    extern float EPSILON_WE,EPSILON_WE_SH;
     extern int EPRECOND_PER_SHOT;
     
     extern int LBFGS_SURFACE;
@@ -509,6 +509,12 @@ void read_par_json(FILE *fp, char *fileinp){
                     }
                     if (get_float_from_objectlist("EPSILON_WE",number_readobjects,&EPSILON_WE,varname_list, value_list))
                         err("Variable EPSILON_WE could not be retrieved from the json input file!");
+                    
+                    if (get_float_from_objectlist("EPSILON_WE_SH",number_readobjects,&EPSILON_WE_SH,varname_list, value_list)) {
+                        EPSILON_WE_SH=EPSILON_WE;
+                        fprintf(fp,"Variable EPSILON_WE_SH is set to EPSILON_WE=%f.\n",EPSILON_WE_SH);
+                    }
+                    
                 }
                 
                 if (get_int_from_objectlist("TESTSHOT_START",number_readobjects,&TESTSHOT_START,varname_list, value_list))