Skip to content

Skip SRF sources with slip_sum=0 and fix sfileoutput #240

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
houjun wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Skip SRF sources with slip_sum=0 and fix sfileoutput #240

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9c7509b
Fix sfile output
houjun Oct 27, 2025
05f926a
Fix sfile output
houjun Oct 27, 2025
a76e0c7
Combine sfile output fix with SRF source updates
houjun Apr 6, 2026
36790fe
Skip creating zero-slip rupture sources
houjun Apr 7, 2026
94f020d
Fix NumPy int alias in HDF5 pytest
houjun Apr 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions src/SfileOutput.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion src/TimeSeries.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
43 changes: 30 additions & 13 deletions src/parseInputFile.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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);
}
Expand Down
43 changes: 30 additions & 13 deletions src/readhdf5.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
Loading