[cig-commits] [commit] master: added displacement gradient output for receivers (e9b2bb0)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Jan 27 15:56:48 PST 2015
Repository : https://github.com/geodynamics/sw4
On branch : master
Link : https://github.com/geodynamics/sw4/compare/d5da240a9d8092ed6aaf96656a305fa284076d53...afc11559354312155d41cd8c298156b3ef67f947
>---------------------------------------------------------------
commit e9b2bb0aa89e1ea22128d09deda826600436983d
Author: Bjorn Sjogreen <sjogreen2 at llnl.gov>
Date: Tue Jan 27 15:53:52 2015 -0800
added displacement gradient output for receivers
>---------------------------------------------------------------
e9b2bb0aa89e1ea22128d09deda826600436983d
src/TimeSeries.C | 93 ++++++++++++++++++++++++++++++++++++++++++++++++++--
src/TimeSeries.h | 2 +-
src/parseInputFile.C | 10 +++++-
src/solve.C | 65 ++++++++++++++++++++++++++++++++++--
4 files changed, 164 insertions(+), 6 deletions(-)
diff --git a/src/TimeSeries.C b/src/TimeSeries.C
index da69a2f..5185a52 100644
--- a/src/TimeSeries.C
+++ b/src/TimeSeries.C
@@ -251,6 +251,8 @@ TimeSeries::TimeSeries( EW* a_ew, std::string fileName, std::string staName, rec
m_nComp=3;
else if (m_mode == Strains)
m_nComp=6;
+ else if (m_mode == DisplacementGradient )
+ m_nComp=9;
// allocate handles to solution array pointers
mRecordedSol = new double*[m_nComp];
@@ -432,7 +434,7 @@ void TimeSeries::writeFile( string suffix )
// get the epicenter from EW object (note that the epicenter is not always known when this object is created)
m_ew->get_epicenter( m_epi_lat, m_epi_lon, m_epi_depth, m_epi_time_offset );
- stringstream ux, uy, uz, uxy, uxz, uyz;
+ stringstream ux, uy, uz, uxy, uxz, uyz, uyx, uzx, uzy;
// Write out displacement components (ux, uy, uz)
@@ -448,7 +450,7 @@ void TimeSeries::writeFile( string suffix )
<< "of size " << mLastTimeStep+1 << ": "
<< filePrefix.str();
- string xfield, yfield, zfield, xyfield, xzfield, yzfield;
+ string xfield, yfield, zfield, xyfield, xzfield, yzfield, yxfield, zxfield, zyfield;
float azimx, azimy, updownang;
if( m_mode == Displacement )
{
@@ -550,6 +552,37 @@ void TimeSeries::writeFile( string suffix )
updownang = 180;
msg << "[xx|yy|zz|xy|xz|yz]" << endl;
}
+ else if( m_mode == DisplacementGradient )
+ {
+ xfield = "DUXDX";
+ xyfield = "DUXDY";
+ xzfield = "DUXDZ";
+
+ yxfield = "DUYDX";
+ yfield = "DUYDY";
+ yzfield = "DUYDZ";
+
+ zxfield = "DUZDX";
+ zyfield = "DUZDY";
+ zfield = "DUZDZ";
+
+ ux << filePrefix.str() << "duxdx";
+ uxy << filePrefix.str() << "duxdy";
+ uxz << filePrefix.str() << "duxdz";
+
+ uyx << filePrefix.str() << "duydx";
+ uy << filePrefix.str() << "duydy";
+ uyz << filePrefix.str() << "duydz";
+
+ uzx << filePrefix.str() << "duzdx";
+ uzy << filePrefix.str() << "duzdy";
+ uz << filePrefix.str() << "duzdz";
+
+ azimx = m_x_azimuth;
+ azimy = m_x_azimuth+90.;
+ updownang = 180;
+ msg << "[duxdx|duxdy|duxdz|duydx|duydy|duydz|duzdx|duzdy|duzdz]" << endl;
+ }
// else if( !m_xycomponent && !m_velocities )
// {
// xfield = "EW";
@@ -729,6 +762,48 @@ void TimeSeries::writeFile( string suffix )
mRecordedFloats[5], (float) m_shift, (float) m_dt,
const_cast<char*>(yzfield.c_str()), 90.0, azimx); // not sure what the updownang or azimuth should be here
}
+ else if (m_mode == DisplacementGradient ) // 9 components
+ {
+ write_sac_format(mLastTimeStep+1,
+ const_cast<char*>(ux.str().c_str()),
+ mRecordedFloats[0], (float) m_shift, (float) m_dt,
+ const_cast<char*>(xfield.c_str()), 90.0, azimx);
+ write_sac_format(mLastTimeStep+1,
+ const_cast<char*>(uxy.str().c_str()),
+ mRecordedFloats[1], (float) m_shift, (float) m_dt,
+ const_cast<char*>(xyfield.c_str()), 90.0, azimx); // not sure what the updownang or azimuth should be here
+ write_sac_format(mLastTimeStep+1,
+ const_cast<char*>(uxz.str().c_str()),
+ mRecordedFloats[2], (float) m_shift, (float) m_dt,
+ const_cast<char*>(xzfield.c_str()), 90.0, azimx); // not sure what the updownang or azimuth should be here
+
+ write_sac_format(mLastTimeStep+1,
+ const_cast<char*>(uyx.str().c_str()),
+ mRecordedFloats[3], (float) m_shift, (float) m_dt,
+ const_cast<char*>(yxfield.c_str()), 90.0, azimx); // not sure what the updownang or azimuth should be here
+ write_sac_format(mLastTimeStep+1,
+ const_cast<char*>(uy.str().c_str()),
+ mRecordedFloats[4], (float) m_shift, (float) m_dt,
+ const_cast<char*>(yfield.c_str()), 90.0, azimy);
+ write_sac_format(mLastTimeStep+1,
+ const_cast<char*>(uyz.str().c_str()),
+ mRecordedFloats[5], (float) m_shift, (float) m_dt,
+ const_cast<char*>(yzfield.c_str()), 90.0, azimx); // not sure what the updownang or azimuth should be here
+
+ write_sac_format(mLastTimeStep+1,
+ const_cast<char*>(uzx.str().c_str()),
+ mRecordedFloats[6], (float) m_shift, (float) m_dt,
+ const_cast<char*>(zxfield.c_str()), 90.0, azimx); // not sure what the updownang or azimuth should be here
+ write_sac_format(mLastTimeStep+1,
+ const_cast<char*>(uzy.str().c_str()),
+ mRecordedFloats[7], (float) m_shift, (float) m_dt,
+ const_cast<char*>(zyfield.c_str()), 90.0, azimx); // not sure what the updownang or azimuth should be here
+ write_sac_format(mLastTimeStep+1,
+ const_cast<char*>(uz.str().c_str()),
+ mRecordedFloats[8], (float) m_shift, (float) m_dt,
+ const_cast<char*>(zfield.c_str()), updownang, 0.0);
+ }
+
} // end if m_sacFormat
if( m_usgsFormat )
@@ -926,6 +1001,18 @@ void TimeSeries::write_usgs_format(string a_fileName)
fprintf(fd, "# Column 6: xz strain component ()\n");
fprintf(fd, "# Column 7: yz strain component ()\n");
}
+ else if( m_mode == DisplacementGradient )
+ {
+ fprintf(fd, "# Column 2 : dux/dx component ()\n");
+ fprintf(fd, "# Column 3 : dux/dy component ()\n");
+ fprintf(fd, "# Column 4 : dux/dz component ()\n");
+ fprintf(fd, "# Column 5 : duy/dx component ()\n");
+ fprintf(fd, "# Column 6 : duy/dy component ()\n");
+ fprintf(fd, "# Column 7 : duy/dz component ()\n");
+ fprintf(fd, "# Column 8 : duz/dx component ()\n");
+ fprintf(fd, "# Column 9 : duz/dy component ()\n");
+ fprintf(fd, "# Column 10: duz/dz component ()\n");
+ }
// write the data
@@ -964,6 +1051,8 @@ void TimeSeries::write_usgs_format(string a_fileName)
printf("strains");
else if( m_mode == Curl )
printf("curl");
+ else if( m_mode == DisplacementGradient )
+ printf("displacement gradient");
printf(" in geographic coordinates\n" );
}
fclose(fd);
diff --git a/src/TimeSeries.h b/src/TimeSeries.h
index 02b6298..22da905 100644
--- a/src/TimeSeries.h
+++ b/src/TimeSeries.h
@@ -47,7 +47,7 @@ class TimeSeries{
public:
// support for derived quantities of the time derivative are not yet implemented
-enum receiverMode{Displacement, Div, Curl, Strains, Velocity /*, DivVelo, CurlVelo, StrainsVelo */ };
+ enum receiverMode{Displacement, Div, Curl, Strains, Velocity, DisplacementGradient /*, DivVelo, CurlVelo, StrainsVelo */ };
TimeSeries( EW* a_ew, std::string fileName, std::string staName, receiverMode mode, bool sacFormat, bool usgsFormat,
double x, double y, double z, bool topoDepth, int writeEvery, bool xyzcomponent=true );
diff --git a/src/parseInputFile.C b/src/parseInputFile.C
index 02e166a..f2d0832 100644
--- a/src/parseInputFile.C
+++ b/src/parseInputFile.C
@@ -4905,7 +4905,8 @@ void EW::processSource(char* buffer, vector<Source*> & a_GlobalUniqueSources )
mxx = myy = mzz = 1;
mxy = mxz = myz = 0;
- bool timereverse = false; // Set true for testing purpose only, the users want to do this themselves outside SW4.
+ bool timereverse = false; // Reverse the SAC data. Set true for testing purpose only, the users want to do this themselves outside SW4.
+ bool useB = false; // Use sac header begin time parameter B.
tDep = iDiscrete6moments;
double dt, t0, latsac, lonsac,cmpazsac, cmpincsac;
@@ -4915,6 +4916,9 @@ void EW::processSource(char* buffer, vector<Source*> & a_GlobalUniqueSources )
fname = basename + ".xx";
bool byteswap;
readSACheader( fname.c_str(), dt, t0, latsac, lonsac, cmpazsac, cmpincsac, utcsac, npts, byteswap );
+ if( !useB )
+ t0 = 0;
+
if( geoCoordSet )
{
double laterr = fabs((latsac-lat)/lat);
@@ -6294,6 +6298,10 @@ void EW::processReceiver(char* buffer, vector<TimeSeries*> & a_GlobalTimeSeries)
{
mode = TimeSeries::Strains;
}
+ else if( strcmp("displacementgradient",token)==0 )
+ {
+ mode = TimeSeries::DisplacementGradient;
+ }
else
{
if (proc_zero())
diff --git a/src/solve.C b/src/solve.C
index 731fa97..d5ac1b5 100644
--- a/src/solve.C
+++ b/src/solve.C
@@ -2751,9 +2751,9 @@ void EW::extractRecordData(TimeSeries::receiverMode mode, int i0, int j0, int k0
// }
}
} // end Curl
- else if(mode == TimeSeries::Strains)
+ else if(mode == TimeSeries::Strains )
{
- uRec.resize(6);
+ uRec.resize(6);
if (g0 < mNumberOfCartesianGrids) // must be a Cartesian grid
{
// int i=m_i0, j=m_j0, k=m_k0, g=m_grid0;
@@ -2806,7 +2806,68 @@ void EW::extractRecordData(TimeSeries::receiverMode mode, int i0, int j0, int k0
uRec[5] = ( 0.5*(duydz+duzdy) );
}
} // end Strains
+ else if(mode == TimeSeries::DisplacementGradient )
+ {
+ uRec.resize(9);
+ if (g0 < mNumberOfCartesianGrids) // must be a Cartesian grid
+ {
+// int i=m_i0, j=m_j0, k=m_k0, g=m_grid0;
+ double factor = 1.0/(2*mGridSize[g0]);
+ double duydx = (U[g0](2,i0+1,j0,k0) - U[g0](2,i0-1,j0,k0))*factor;
+ double duzdx = (U[g0](3,i0+1,j0,k0) - U[g0](3,i0-1,j0,k0))*factor;
+ double duxdy = (U[g0](1,i0,j0+1,k0) - U[g0](1,i0,j0-1,k0))*factor;
+ double duzdy = (U[g0](3,i0,j0+1,k0) - U[g0](3,i0,j0-1,k0))*factor;
+ double duxdz = (U[g0](1,i0,j0,k0+1) - U[g0](1,i0,j0,k0-1))*factor;
+ double duydz = (U[g0](2,i0,j0,k0+1) - U[g0](2,i0,j0,k0-1))*factor;
+ double duxdx = (U[g0](1,i0+1,j0,k0) - U[g0](1,i0-1,j0,k0))*factor;
+ double duydy = (U[g0](2,i0,j0+1,k0) - U[g0](2,i0,j0-1,k0))*factor;
+ double duzdz = (U[g0](3,i0,j0,k0+1) - U[g0](3,i0,j0,k0-1))*factor;
+ uRec[0] = duxdx;
+ uRec[1] = duxdy;
+ uRec[2] = duxdz;
+ uRec[3] = duydx;
+ uRec[4] = duydy;
+ uRec[5] = duydz;
+ uRec[6] = duzdx;
+ uRec[7] = duzdy;
+ uRec[8] = duzdz;
+ }
+ else // must be curvilinear
+ {
+// int i=m_i0, j=m_j0, k0=m_k00, g0=m_grid0;
+ double factor = 0.5/sqrt(mJ(i0,j0,k0));
+ double duxdq = (U[g0](1,i0+1,j0,k0) - U[g0](1,i0-1,j0,k0));
+ double duydq = (U[g0](2,i0+1,j0,k0) - U[g0](2,i0-1,j0,k0));
+ double duzdq = (U[g0](3,i0+1,j0,k0) - U[g0](3,i0-1,j0,k0));
+ double duxdr = (U[g0](1,i0,j0+1,k0) - U[g0](1,i0,j0-1,k0));
+ double duydr = (U[g0](2,i0,j0+1,k0) - U[g0](2,i0,j0-1,k0));
+ double duzdr = (U[g0](3,i0,j0+1,k0) - U[g0](3,i0,j0-1,k0));
+ double duxds = (U[g0](1,i0,j0,k0+1) - U[g0](1,i0,j0,k0-1));
+ double duyds = (U[g0](2,i0,j0,k0+1) - U[g0](2,i0,j0,k0-1));
+ double duzds = (U[g0](3,i0,j0,k0+1) - U[g0](3,i0,j0,k0-1));
+ double duzdy = (mMetric(1,i0,j0,k0)*duzdr+mMetric(3,i0,j0,k0)*duzds)*factor;
+ double duydz = (mMetric(4,i0,j0,k0)*duyds)*factor;
+ double duxdz = (mMetric(4,i0,j0,k0)*duxds)*factor;
+ double duzdx = (mMetric(1,i0,j0,k0)*duzdq+mMetric(2,i0,j0,k0)*duzds)*factor;
+ double duydx = (mMetric(1,i0,j0,k0)*duydq+mMetric(2,i0,j0,k0)*duyds)*factor;
+ double duxdy = (mMetric(1,i0,j0,k0)*duxdr+mMetric(3,i0,j0,k0)*duxds)*factor;
+ double duxdx = ( mMetric(1,i0,j0,k0)*(U[g0](1,i0+1,j0,k0) - U[g0](1,i0-1,j0,k0))+
+ mMetric(2,i0,j0,k0)*(U[g0](1,i0,j0,k0+1) - U[g0](1,i0,j0,k0-1)) )*factor;
+ double duydy = ( mMetric(1,i0,j0,k0)*(U[g0](2,i0,j0+1,k0) - U[g0](2,i0,j0-1,k0))+
+ mMetric(3,i0,j0,k0)*(U[g0](2,i0,j0,k0+1) - U[g0](2,i0,j0,k0-1)) )*factor;
+ double duzdz = ( mMetric(4,i0,j0,k0)*(U[g0](3,i0,j0,k0+1) - U[g0](3,i0,j0,k0-1)) )*factor;
+ uRec[0] = duxdx;
+ uRec[1] = duxdy;
+ uRec[2] = duxdz;
+ uRec[3] = duydx;
+ uRec[4] = duydy;
+ uRec[5] = duydz;
+ uRec[6] = duzdx;
+ uRec[7] = duzdy;
+ uRec[8] = duzdz;
+ }
+ } // end DisplacementGradient
return;
}
More information about the CIG-COMMITS
mailing list