diff --git a/pytest/verify_hdf5.py b/pytest/verify_hdf5.py index 703d1106..d7690a47 100644 --- a/pytest/verify_hdf5.py +++ b/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] diff --git a/src/SfileOutput.C b/src/SfileOutput.C index ad28fa7a..0344e1c8 100644 --- a/src/SfileOutput.C +++ b/src/SfileOutput.C @@ -739,6 +739,14 @@ void SfileOutput::write_image(const char *fname, std::vector& 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); diff --git a/src/TimeSeries.C b/src/TimeSeries.C index eea03c6f..3be13659 100644 --- a/src/TimeSeries.C +++ b/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) diff --git a/src/parseInputFile.C b/src/parseInputFile.C index 21d2a31b..919e8208 100644 --- a/src/parseInputFile.C +++ b/src/parseInputFile.C @@ -6263,7 +6263,7 @@ void EW::processRupture(char* buffer, vector > & 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 > & 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 > & 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 > & 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); } diff --git a/src/readhdf5.C b/src/readhdf5.C index fd728b52..85151861 100644 --- a/src/readhdf5.C +++ b/src/readhdf5.C @@ -569,7 +569,7 @@ void readRuptureHDF5(char *fname, vector > & 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 > & 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 > & 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 > & 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)