Skip SRF sources with slip_sum=0 and fix sfileoutput

pytest/verify_hdf5.py

@@ -79,10 +79,10 @@ def read_sw4img(fname):
     # Grid size info
     h = np.zeros(npatch, dtype=np.float64)
     zmin = np.zeros(npatch, dtype=np.float64)
-    i = np.zeros(npatch, dtype=np.int)
-    ni = np.zeros(npatch, dtype=np.int)
-    j = np.zeros(npatch, dtype=np.int)
-    nj = np.zeros(npatch, dtype=np.int)
+    i = np.zeros(npatch, dtype=int)
+    ni = np.zeros(npatch, dtype=int)
+    j = np.zeros(npatch, dtype=int)
+    nj = np.zeros(npatch, dtype=int)
     nelem = 0
     for u in range (0, npatch):
         h[u] = struct.unpack('d', img.read(8))[0]

src/SfileOutput.C

@@ -739,6 +739,14 @@ void SfileOutput::write_image(const char *fname, std::vector<Sarray>& a_Z )
     H5Awrite(attr, H5T_NATIVE_DOUBLE, &spacing);
     H5Aclose(attr);
 
+    aname = "Coarsest horizontal grid spacing";
+    spacing = mEW->mGridSize[0]*stH;
+    attr = H5Acreate(h5_fid, aname, H5T_NATIVE_DOUBLE, attr_space1, H5P_DEFAULT, H5P_DEFAULT);
+    if( attr < 0 )
+      VERIFY2(0, "ERROR: SfileOutput::write_image, error creating " << aname);
+    H5Awrite(attr, H5T_NATIVE_DOUBLE, &spacing);
+    H5Aclose(attr);
+
     aname = "Attenuation";
     int att = mEW->usingAttenuation();
     attr = H5Acreate(h5_fid, aname, H5T_NATIVE_INT, attr_space1, H5P_DEFAULT, H5P_DEFAULT);

src/TimeSeries.C

@@ -1329,7 +1329,7 @@ write_hdf5_format(int npts, hid_t grp, float *y, float btime, float dt, char *va
 
   if (isLast && ret == 1) {
     m_nptsWritten += count;
-    H5Gflush(grp);
+    // H5Gflush(grp);
   }
 
   if (mDownSample > 1) 

src/parseInputFile.C

@@ -6263,7 +6263,7 @@ void EW::processRupture(char* buffer, vector<vector<Source*> > & a_GlobalUniqueS
       printf("Number of point sources in data block: %i\n", npts);
 
 // read all point sources
-    int nSources=0, nu1=0, nu2=0, nu3=0;
+    int nSources=0, nu1=0, nu2=0, nu3=0, nskip_zero_slip=0;
     for (int pts=0; pts<npts; pts++) 
     {
       double lon, lat, dep, stk, dip, area, tinit, dt, rake, slip1, slip2, slip3;
@@ -6342,20 +6342,35 @@ void EW::processRupture(char* buffer, vector<vector<Source*> > & a_GlobalUniqueS
         {
            printf("INFO: SRF file: dt*sum(slip_vel)=%e [m], total slip (from header)=%e [m]\n", slip_sum, slip_m);
         }
-// scale time series to sum to integrate to one        
-	for (int i=1; i<=nt1dim+1; i++)
-	{
-           par[i] /= slip_sum;
-	}
-        if (proc_zero() && mVerbose >= 2)
+        float_sw4 slip_sum_tol = 1e-12;
+        bool skip_zero_slip_point = false;
+        if( slip_sum > -slip_sum_tol && slip_sum < slip_sum_tol )
         {
-           slip_sum=0;
-           for (int i=1; i<=nt1dim+1; i++)
+           nskip_zero_slip++;
+           skip_zero_slip_point = true;
+           if( proc_zero() && nskip_zero_slip <= 10 )
            {
-              slip_sum += par[i];
+              printf("WARNING: skipping rupture point #%i because dt*sum(slip_vel)=%e [m], slip1=%e [m]\n",
+                     pts+1, slip_sum, slip_m);
            }
-           slip_sum *=dt;
-           printf("INFO: SRF file: After scaling time series: dt*sum(par)=%e [m]\n", slip_sum);
+        }
+// scale time series to sum to integrate to one        
+	if( !skip_zero_slip_point )
+	{
+	  for (int i=1; i<=nt1dim+1; i++)
+	  {
+             par[i] /= slip_sum;
+	  }
+          if (proc_zero() && mVerbose >= 2)
+          {
+             slip_sum=0;
+             for (int i=1; i<=nt1dim+1; i++)
+             {
+                slip_sum += par[i];
+             }
+             slip_sum *=dt;
+             printf("INFO: SRF file: After scaling time series: dt*sum(par)=%e [m]\n", slip_sum);
+          }
         }
 //done scaling        
 
@@ -6433,7 +6448,7 @@ void EW::processRupture(char* buffer, vector<vector<Source*> > & a_GlobalUniqueS
 	  if (m_myRank == 0)
 	    cout << sourceposerr.str();
 	}
-	else if( event_is_in_proc(event) )
+	else if( !skip_zero_slip_point && event_is_in_proc(event) )
 	{
            event = global_to_local_event(event);
            sourcePtr = new Source(this, freq, t0, x, y, z, mxx, mxy, mxz, myy, myz, mzz,
@@ -6511,6 +6526,8 @@ void EW::processRupture(char* buffer, vector<vector<Source*> > & a_GlobalUniqueS
     if (proc_zero())
       printf("Read npts=%i, made %i point moment tensor sources, nu1=%i, nu2=%i, nu3=%i\n", 
 	     npts, nSources, nu1, nu2, nu3);
+    if (proc_zero() && nskip_zero_slip > 0)
+      printf("Skipped %i rupture points with zero slip-velocity integral in u1.\n", nskip_zero_slip);
 
     fclose(fd);
   }

src/readhdf5.C

@@ -569,7 +569,7 @@ void readRuptureHDF5(char *fname, vector<vector<Source*> > & a_GlobalUniqueSourc
   int npts = 0, nseg = 0, nsr1 = 0;
   hsize_t dims;
   double rVersion;
-  int nSources=0, nu1=0, nu2=0, nu3=0;
+  int nSources=0, nu1=0, nu2=0, nu3=0, nskip_zero_slip=0;
 
   stime = MPI_Wtime();
   // Only rank 0 reads data, then broadcast to all other processes
@@ -777,20 +777,35 @@ void readRuptureHDF5(char *fname, vector<vector<Source*> > & a_GlobalUniqueSourc
       {
          printf("INFO: SRF file: dt*sum(slip_vel)=%e [m], total slip (from header)=%e [m]\n", slip_sum, slip_m);
       }
-      // scale time series to sum to integrate to one        
-      for (int i=1; i<=nt1dim+1; i++)
+      float_sw4 slip_sum_tol = 1e-12;
+      bool skip_zero_slip_point = false;
+      if( slip_sum > -slip_sum_tol && slip_sum < slip_sum_tol )
       {
-         par[i] /= slip_sum;
+        nskip_zero_slip++;
+        skip_zero_slip_point = true;
+        if( world_rank == 0 && nskip_zero_slip <= 10 )
+        {
+          printf("WARNING: skipping rupture point #%i because dt*sum(slip_vel)=%e [m], slip1=%e [m]\n",
+                 pts+1, slip_sum, slip_m);
+        }
       }
-      if (world_rank == 0 && mVerbose >= 2)
+      // scale time series to sum to integrate to one
+      if( !skip_zero_slip_point )
       {
-         slip_sum=0;
-         for (int i=1; i<=nt1dim+1; i++)
-         {
-            slip_sum += par[i];
-         }
-         slip_sum *=dt;
-         printf("INFO: SRF file: After scaling time series: dt*sum(par)=%e [m]\n", slip_sum);
+        for (int i=1; i<=nt1dim+1; i++)
+        {
+           par[i] /= slip_sum;
+        }
+        if (world_rank == 0 && mVerbose >= 2)
+        {
+           slip_sum=0;
+           for (int i=1; i<=nt1dim+1; i++)
+           {
+              slip_sum += par[i];
+           }
+           slip_sum *=dt;
+           printf("INFO: SRF file: After scaling time series: dt*sum(par)=%e [m]\n", slip_sum);
+        }
       }
       //done scaling        
 
@@ -868,7 +883,7 @@ void readRuptureHDF5(char *fname, vector<vector<Source*> > & a_GlobalUniqueSourc
         if (world_rank == 0)
           cout << sourceposerr.str();
       }
-      else
+      else if( !skip_zero_slip_point )
       {
         sourcePtr = new Source(ew, freq, t0, x, y, z, mxx, mxy, mxz, myy, myz, mzz,
                                tDep, formstring, topodepth, ncyc, par, npar, ipar, nipar, true ); // true is correctStrengthForMu
@@ -911,6 +926,8 @@ void readRuptureHDF5(char *fname, vector<vector<Source*> > & a_GlobalUniqueSourc
   } // end for all sources
   if (world_rank == 0)
     printf("Read npts=%i, made %i point moment tensor sources, nu1=%i, nu2=%i, nu3=%i\n", npts, nSources, nu1, nu2, nu3);
+  if (world_rank == 0 && nskip_zero_slip > 0)
+    printf("Skipped %i rupture points with zero slip-velocity integral in u1.\n", nskip_zero_slip);
 
   etime = MPI_Wtime();
   if (is_debug && world_rank == 0)