[cig-commits] r14264 - short/2.5D/benchmarks/savageprescott/utils
willic3 at geodynamics.org
willic3 at geodynamics.org
Sun Mar 8 19:44:28 PDT 2009
Author: willic3
Date: 2009-03-08 19:44:28 -0700 (Sun, 08 Mar 2009)
New Revision: 14264
Modified:
short/2.5D/benchmarks/savageprescott/utils/README
short/2.5D/benchmarks/savageprescott/utils/savpres_ss.cfg
short/2.5D/benchmarks/savageprescott/utils/savpres_ss.py
Log:
Updated analytical solution to compute velocity results, and updated
README to explain about generating and comparing velocity results.
Modified: short/2.5D/benchmarks/savageprescott/utils/README
===================================================================
--- short/2.5D/benchmarks/savageprescott/utils/README 2009-03-09 02:11:03 UTC (rev 14263)
+++ short/2.5D/benchmarks/savageprescott/utils/README 2009-03-09 02:44:28 UTC (rev 14264)
@@ -1,12 +1,24 @@
-This directory contains a Python code to compute the analytical solution
-for the Savage and Prescott [1978] model, along with a .cfg file to specify
-the problem. To run the model, simply type './savpres_ss.py'. The
-savpres_ss.cfg is read automatically. Once finished, vtk results will be
-stored in the 'vtk_results' directory, and a CSV file will also be produced
-(savageprescott.csv). The CSV file is meant to be used with a plotting
-package. The vtk results probably are not very useful.
+This directory contains two Python utility codes, along with the .cfg files
+needed to run them. The first code (savpres_ss.py) computes the analytical
+solution for displacement and velocity for the Savage and Prescott [1978]
+model. To run the model, simply type:
-Note that to compare the results against those from PyLith, you should
-reference the displacements to the fist step in a cycle. For example, to
-look at the results from cycle 5, you would subtract the results in column
-c5_t0000 from all of the other columns in cycle 5.
+./savpres_ss.py
+
+The savpres_ss.cfg file is read automatically. At present, things are set
+up only for CSV (comma-separated values) output, since the VTK output is
+not likely to be useful. The CSV results can be used with most plotting
+packages.
+
+Note that you will not be able to directly compare the analytical and
+numerical solutions (see top-level README). You can instead compare the
+velocity results for time steps in the later cycles.
+
+The second utility code (vtkdiff.py) allows you to compute velocities from
+the displacement results produced by PyLith. To do this, you just type:
+
+./vtkdiff.py vtkdiff_hex8_unif1_20km.cfg
+
+This will generate velocity results for the uniform mesh with 20 km
+resolution. Results are placed in the results directory for the
+appropriate model.
Modified: short/2.5D/benchmarks/savageprescott/utils/savpres_ss.cfg
===================================================================
--- short/2.5D/benchmarks/savageprescott/utils/savpres_ss.cfg 2009-03-09 02:11:03 UTC (rev 14263)
+++ short/2.5D/benchmarks/savageprescott/utils/savpres_ss.cfg 2009-03-09 02:44:28 UTC (rev 14264)
@@ -9,13 +9,22 @@
shear_modulus = 30.0*GPa
plate_velocity = 2.0*cm/year
number_cycles = 10
-number_steps = 20
+number_steps = 40
number_terms = 100
number_points = 400
delta_x = 5.0*km
x_epsilon = 0.01*m
-vtk_basename = vtk_results/savageprescott.vtk
-csv_filename = savageprescott.csv
+output_displ_vtk = False
+output_vel_vtk = False
+output_displ_csv = True
+output_vel_csv = True
+# displ_vtk_basename = vtk_results/savpres_displ.vtk
+displ_csv_filename = savpres_displ.csv
+# vel_vtk_basename = vtk_results/savpres_vel.vtk
+vel_csv_filename = savpres_vel.csv
+# Convert meters to cm and m/s to cm/year
+displ_scale_factor = 100.0
+vel_scale_factor = 3.15576e9
time_units = 1.0*year
time_stamp_width = 4
title = Test model with D = 20 km, H = 40 km, and tau0 = 0.5.
Modified: short/2.5D/benchmarks/savageprescott/utils/savpres_ss.py
===================================================================
--- short/2.5D/benchmarks/savageprescott/utils/savpres_ss.py 2009-03-09 02:11:03 UTC (rev 14263)
+++ short/2.5D/benchmarks/savageprescott/utils/savpres_ss.py 2009-03-09 02:44:28 UTC (rev 14264)
@@ -49,8 +49,16 @@
## @li \b number_points Number of points at which to compute solution.
## @li \b delta_x Distance between computation points.
## @li \b x_epsilon Offset for computation point closest to the fault.
- ## @li \b vtk_basename Base name for VTK output files.
- ## @li \b csv_filename Filename for CSV output.
+ ## @li \b output_displ_vtk Output displacement VTK results?
+ ## @li \b output_displ_csv Output displacement CSV results?
+ ## @li \b output_vel_vtk Output velocity VTK results?
+ ## @li \b output_vel_csv Output velocity CSV results?
+ ## @li \b displ_vtk_basename Base name for VTK displacement output files.
+ ## @li \b displ_csv_filename Filename for CSV displacement output.
+ ## @li \b vel_vtk_basename Base name for VTK velocity output files.
+ ## @li \b vel_csv_filename Filename for CSV velocity output.
+ ## @li \b displ_scale_factor Scaling factor for output displacements.
+ ## @li \b vel_scale_factor Scaling factor for output velocities.
## @li \b time_units Time units to use for output filenames.
## @li \b time_stamp_width Width of time stamp in output filenames.
## @li \b title Title to appear at the top of VTK files.
@@ -103,12 +111,40 @@
xEpsilon = pyre.inventory.dimensional("x_epsilon", default=0.001*m)
xEpsilon.meta['tip'] = "Offset for computation point closest to the fault."
- VTKBaseName = pyre.inventory.str("vtk_basename", default="savpres_ss.vtk")
- VTKBaseName.meta['tip'] = "Base filename of VTK output."
+ outputDisplVTK = pyre.inventory.bool("output_displ_vtk", default=False)
+ outputDisplVTK.meta['tip'] = "Output displacement VTK files?"
- CSVFileName = pyre.inventory.str("csv_filename", default="savpres_ss.csv")
- CSVFileName.meta['tip'] = "Filename for CSV output."
+ outputDisplCSV = pyre.inventory.bool("output_displ_csv", default=True)
+ outputDisplCSV.meta['tip'] = "Output displacement CSV file?"
+ outputVelVTK = pyre.inventory.bool("output_vel_vtk", default=False)
+ outputVelVTK.meta['tip'] = "Output velocity VTK files?"
+
+ outputVelCSV = pyre.inventory.bool("output_vel_csv", default=True)
+ outputVelCSV.meta['tip'] = "Output velocity CSV file?"
+
+ displVTKBaseName = pyre.inventory.str("displ_vtk_basename",
+ default="savpres_ss_displ.vtk")
+ displVTKBaseName.meta['tip'] = "Base filename of VTK displacement output."
+
+ displCSVFileName = pyre.inventory.str("displ_csv_filename",
+ default="savpres_ss_displ.csv")
+ displCSVFileName.meta['tip'] = "Filename for CSV displacement output."
+
+ velVTKBaseName = pyre.inventory.str("vel_vtk_basename",
+ default="savpres_ss_vel.vtk")
+ velVTKBaseName.meta['tip'] = "Base filename of VTK velocity output."
+
+ velCSVFileName = pyre.inventory.str("vel_csv_filename",
+ default="savpres_ss_vel.csv")
+ velCSVFileName.meta['tip'] = "Filename for CSV velocity output."
+
+ displScaleFactor = pyre.inventory.float("displ_scale_factor", default=1.0)
+ displScaleFactor.meta['tip'] = "Scale factor for displacement output."
+
+ velScaleFactor = pyre.inventory.float("vel_scale_factor", default=1.0)
+ velScaleFactor.meta['tip'] = "Scale factor for velocity output."
+
timeUnits = pyre.inventory.dimensional("time_units", default=1.0*year)
timeUnits.meta['tip'] = "Time units to use for output filenames."
@@ -131,8 +167,14 @@
# pdb.set_trace()
self._genPoints()
self._genSolution()
- self._writeSolutionVTK()
- self._writeSolutionCSV()
+ if self.outputDisplVTK:
+ self._writeSolutionVTK("displacement")
+ if self.outputVelVTK:
+ self._writeSolutionVTK("velocity")
+ if self.outputDisplCSV:
+ self._writeSolutionCSV("displacement")
+ if self.outputVelCSV:
+ self._writeSolutionCSV("velocity")
return
@@ -155,8 +197,16 @@
self.numberPoints = self.inventory.numberPoints
self.deltaX = self.inventory.deltaX.value
self.xEpsilon = self.inventory.xEpsilon.value
- self.VTKBaseName = self.inventory.VTKBaseName
- self.CSVFileName = self.inventory.CSVFileName
+ self.outputDisplVTK = self.inventory.outputDisplVTK
+ self.outputDisplCSV = self.inventory.outputDisplCSV
+ self.outputVelVTK = self.inventory.outputVelVTK
+ self.outputVelCSV = self.inventory.outputVelCSV
+ self.displVTKBaseName = self.inventory.displVTKBaseName
+ self.displCSVFileName = self.inventory.displCSVFileName
+ self.velVTKBaseName = self.inventory.velVTKBaseName
+ self.velCSVFileName = self.inventory.velCSVFileName
+ self.displScaleFactor = self.inventory.displScaleFactor
+ self.velScaleFactor = self.inventory.velScaleFactor
self.timeUnits = self.inventory.timeUnits.value
self.timeStampWidth = self.inventory.timeStampWidth
self.title = self.inventory.title
@@ -173,9 +223,9 @@
Create array of points for output along with series terms
for each point.
"""
- self.points = numpy.zeros(self.numberPoints, dtype=float)
+ self.points = numpy.zeros(self.numberPoints, dtype=numpy.float64)
self.pointCoeff = numpy.zeros((self.numberPoints, self.numberTerms),
- dtype=float)
+ dtype=numpy.float64)
for point in range(self.numberPoints):
self.points[point] = max(self.xEpsilon, point*self.deltaX)
@@ -200,32 +250,50 @@
solutionU2 = numpy.zeros((self.numberCycles,
self.numberSteps + 1,
self.numberPoints),
- dtype=float)
+ dtype=numpy.float64)
self.solutionUTot = numpy.zeros((self.numberCycles,
self.numberSteps + 1,
self.numberPoints),
- dtype=float)
- oneArray = numpy.ones(self.numberPoints, dtype=float)
+ dtype=numpy.float64)
+ solutionV2 = numpy.zeros((self.numberCycles,
+ self.numberSteps + 1,
+ self.numberPoints),
+ dtype=numpy.float64)
+ self.solutionVTot = numpy.zeros((self.numberCycles,
+ self.numberSteps + 1,
+ self.numberPoints),
+ dtype=numpy.float64)
+ oneArray = numpy.ones(self.numberPoints, dtype=numpy.float64)
for cycle in range(self.numberCycles):
time = cycle * self.numberSteps * self.deltaT
tau = time * self.tauFac
if cycle > 0:
solutionU2[cycle, :, :] += solutionU2[cycle - 1, :, :]
+ solutionV2[cycle, :, :] += solutionV2[cycle - 1, :, :]
for step in range(self.numberSteps + 1):
if cycle == 0:
- solutionT = self._u2A(tau)
+ solutionUT = self._u2A(tau)[0]
+ solutionVT = self._u2A(tau)[1]
else:
- solutionT = self._u2B(tau)
+ solutionUT = self._u2B(tau)[0]
+ solutionVT = self._u2B(tau)[1]
- solutionU2[cycle, step, :] += solutionT
+ solutionU2[cycle, step, :] += solutionUT
+ solutionV2[cycle, step, :] += solutionVT
self.solutionUTot[cycle, step, :] = solutionU2[cycle, step, :] + \
time * self.velocity * oneArray
+ self.solutionVTot[cycle, step, :] = self.tauFac * \
+ solutionV2[cycle, step, :] + \
+ self.velocity * oneArray
time = time + self.deltaT
tau = time * self.tauFac
+ self.solutionUTot *= self.displScaleFactor
+ self.solutionVTot *= self.velScaleFactor
+
return
@@ -249,31 +317,36 @@
"""
Computes viscoelastic solution for times less than the recurrence time.
"""
- solution = numpy.zeros(self.numberPoints, dtype=float)
+ solutionU = numpy.zeros(self.numberPoints, dtype=numpy.float64)
+ solutionV = numpy.zeros(self.numberPoints, dtype=numpy.float64)
for point in range(self.numberPoints):
- solution[point] = tau*(-0.5 * math.pi + \
- numpy.arctan(self.points[point]/self.lockDepth))/ \
- self.tau0
+ solution = (-0.5 * math.pi + \
+ numpy.arctan(self.points[point]/self.lockDepth))/self.tau0
+ solutionU[point] = tau * solution
+ solutionV[point] = solution
aPrev = 0.0
bPrev = 0.0
factPrev = 1.0
for term in range(self.numberTerms):
aN, bN, factN = self._timeCoeff(term, tau, aPrev, bPrev, factPrev)
- solution[point] -= bN * self.pointCoeff[point, term]
+ solutionU[point] -= bN * self.pointCoeff[point, term]
+ solutionV[point] -= aN * self.pointCoeff[point, term]/self.tau0
aPrev = aN
bPrev = bN
factPrev = factN
- solution *= 2.0 * self.velocity * self.recurrenceTime/math.pi
- return solution
+ solutionU *= 2.0 * self.velocity * self.recurrenceTime/math.pi
+ solutionV *= 2.0 * self.velocity * self.recurrenceTime/math.pi
+ return [solutionU, solutionV]
def _u2B(self, tau):
"""
Computes viscoelastic solution for times greater than the recurrence time.
"""
- solution = numpy.zeros(self.numberPoints, dtype=float)
+ solutionU = numpy.zeros(self.numberPoints, dtype=numpy.float64)
+ solutionV = numpy.zeros(self.numberPoints, dtype=numpy.float64)
tau2 = tau - self.tau0
for point in range(self.numberPoints):
@@ -287,8 +360,11 @@
a1N, b1N, fact1N = self._timeCoeff(term, tau, a1Prev, b1Prev, fact1Prev)
a2N, b2N, fact2N = self._timeCoeff(term, tau2, a2Prev, b2Prev,
fact2Prev)
- solution[point] += self.pointCoeff[point, term] * \
+ daDt = tau2**term * math.exp(-tau2)/fact2N
+ solutionU[point] += self.pointCoeff[point, term] * \
(b2N - b1N + a2N)
+ solutionV[point] += self.pointCoeff[point, term] * \
+ (a2N/self.tau0 - a1N/self.tau0 + daDt)
a1Prev = a1N
b1Prev = b1N
fact1Prev = fact1N
@@ -296,22 +372,30 @@
b2Prev = b2N
fact2Prev = fact2N
- solution *= 2.0 * self.velocity * self.recurrenceTime/math.pi
- return solution
+ solutionU *= 2.0 * self.velocity * self.recurrenceTime/math.pi
+ solutionV *= 2.0 * self.velocity * self.recurrenceTime/math.pi
+ return [solutionU, solutionV]
- def _writeSolutionVTK(self):
+ def _writeSolutionVTK(self, solutionType):
"""
Generate VTK filename and write results to file.
"""
- if self.VTKBaseName.endswith('.vtk'):
- fileBase = self.VTKBaseName[:self.VTKBaseName.rfind('.vtk')]
- elif self.VTKBaseName.endswith('.VTK'):
- fileBase = self.VTKBaseName[:self.VTKBaseName.rfind('.VTK')]
+ if solutionType == "displacement":
+ VTKBaseName = self.displVTKBaseName
+ solution = self.solutionUTot
else:
- fileBase = self.VTKBaseName
+ VTKBaseName = self.velVTKBaseName
+ solution = self.solutionVTot
+ if VTKBaseName.endswith('.vtk'):
+ fileBase = VTKBaseName[:VTKBaseName.rfind('.vtk')]
+ elif VTKBaseName.endswith('.VTK'):
+ fileBase = VTKBaseName[:VTKBaseName.rfind('.VTK')]
+ else:
+ fileBase = VTKBaseName
+
for cycle in range(self.numberCycles):
fileBaseCycle = fileBase + "_c" + str(cycle) + "_t"
time = 0.0
@@ -321,14 +405,14 @@
timeStampString = repr(timeStampInt).rjust(self.timeStampWidth, '0')
VTKFile = fileBaseCycle + timeStampString + ".vtk"
f = open(VTKFile, 'w')
- self._writeVTK(f, cycle, step)
+ self._writeVTK(f, solution, solutionType, cycle, step)
f.close()
time += self.deltaT
return
- def _writeVTK(self, f, cycle, step):
+ def _writeVTK(self, f, solution, solutionType, cycle, step):
"""
Writes solution to VTK file as a set of points.
"""
@@ -344,23 +428,30 @@
f.write('\n')
f.write('POINT_DATA '+str(self.numberPoints)+'\n')
- f.write('SCALARS displacement double 3\n')
+ f.write('SCALARS '+solutionType+' double 3\n')
f.write('LOOKUP_TABLE default\n')
uX = 0.0
uZ = 0.0
for point in range(self.numberPoints):
f.write(' %.12g %.12g %.12g\n' %
- (uX, self.solutionUTot[cycle, step, point], uZ))
+ (uX, solution[cycle, step, point], uZ))
return
- def _writeSolutionCSV(self):
+ def _writeSolutionCSV(self, solutionType):
"""
Write solution to a CSV file.
"""
- f = open(self.CSVFileName, 'w')
+ if solutionType == "displacement":
+ CSVFileName = self.displCSVFileName
+ solution = self.solutionUTot
+ else:
+ CSVFileName = self.velCSVFileName
+ solution = self.solutionVTot
+
+ f = open(CSVFileName, 'w')
head = "X,Y,Z"
for cycle in range(self.numberCycles):
cycleHead = "c" + str(cycle) + "_t"
@@ -380,7 +471,7 @@
f.write(' %.12g, %.12g, %.12g' % (self.points[point], y, z))
for cycle in range(self.numberCycles):
for step in range(self.numberSteps + 1):
- f.write(', %.12g' % self.solutionUTot[cycle, step, point])
+ f.write(', %.12g' % solution[cycle, step, point])
f.write('\n')
f.close()
More information about the CIG-COMMITS
mailing list