[cig-commits] r14326 - short/3D/PyLith/trunk/playpen/savpres_ss
willic3 at geodynamics.org
willic3 at geodynamics.org
Sun Mar 15 17:21:01 PDT 2009
Author: willic3
Date: 2009-03-15 17:21:01 -0700 (Sun, 15 Mar 2009)
New Revision: 14326
Modified:
short/3D/PyLith/trunk/playpen/savpres_ss/savpres_ss.cfg
short/3D/PyLith/trunk/playpen/savpres_ss/savpres_ss.py
Log:
Included ability to compute velocities as well as displacements, and added
many more options.
These versions are now identical to the versions in the
2.5D/benchmarks/savageprescott/utils directory.
Modified: short/3D/PyLith/trunk/playpen/savpres_ss/savpres_ss.cfg
===================================================================
--- short/3D/PyLith/trunk/playpen/savpres_ss/savpres_ss.cfg 2009-03-16 00:16:26 UTC (rev 14325)
+++ short/3D/PyLith/trunk/playpen/savpres_ss/savpres_ss.cfg 2009-03-16 00:21:01 UTC (rev 14326)
@@ -5,17 +5,34 @@
elas_thick = 40.0*km
lock_depth = 20.0*km
recurrence_time = 200.0*year
-viscosity = 9.46728e19*Pa*s
+viscosity = 4.73364e+19*Pa*s
shear_modulus = 30.0*GPa
-plate_velocity = 1.0*cm/year
+plate_velocity = 2.0*cm/year
+
number_cycles = 10
-number_steps = 10
-number_terms = 50
-number_points = 100
+number_steps = 40
+number_terms = 100
+number_points = 400
delta_x = 5.0*km
x_epsilon = 0.01*m
-vtk_basename = 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
+
+# Use km for output coordinates
+coord_scale_factor = 0.001
+coord_units = km
+
time_units = 1.0*year
time_stamp_width = 4
-title = Test model with D = 20 km, H = 40 km, and tau0 = 1.0.
+title = Test model with D = 20 km, H = 40 km, and tau0 = 0.5.
Modified: short/3D/PyLith/trunk/playpen/savpres_ss/savpres_ss.py
===================================================================
--- short/3D/PyLith/trunk/playpen/savpres_ss/savpres_ss.py 2009-03-16 00:16:26 UTC (rev 14325)
+++ short/3D/PyLith/trunk/playpen/savpres_ss/savpres_ss.py 2009-03-16 00:21:01 UTC (rev 14326)
@@ -49,8 +49,18 @@
## @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 coord_scale_factor Scaling factor for output coordinates.
+ ## @li \b vel_scale_factor Scaling factor for output velocities.
+ ## @li \b coord_units Units used for output coordinates.
## @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 +113,46 @@
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."
+
+ coordScaleFactor = pyre.inventory.float("coord_scale_factor", default=1.0)
+ coordScaleFactor.meta['tip'] = "Scale factor for output coordinates."
+
+ coordUnits = pyre.inventory.str("coord_units", default="m")
+ coordUnits.meta['tip'] = "Units used for output coordinates."
+
timeUnits = pyre.inventory.dimensional("time_units", default=1.0*year)
timeUnits.meta['tip'] = "Time units to use for output filenames."
@@ -131,8 +175,15 @@
# pdb.set_trace()
self._genPoints()
self._genSolution()
- self._writeSolutionVTK()
- self._writeSolutionCSV()
+ self.points *= self.coordScaleFactor
+ 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 +206,18 @@
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.coordScaleFactor = self.inventory.coordScaleFactor
+ self.coordUnits = self.inventory.coordUnits
self.timeUnits = self.inventory.timeUnits.value
self.timeStampWidth = self.inventory.timeStampWidth
self.title = self.inventory.title
@@ -173,9 +234,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 +261,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 +328,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 +371,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 +383,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 +416,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,43 +439,48 @@
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')
- head = "X,Y,Z"
+ if solutionType == "displacement":
+ CSVFileName = self.displCSVFileName
+ solution = self.solutionUTot
+ else:
+ CSVFileName = self.velCSVFileName
+ solution = self.solutionVTot
+
+ f = open(CSVFileName, 'w')
+ head = "Distance from Fault (" + self.coordUnits + ")"
for cycle in range(self.numberCycles):
- cycleHead = "c" + str(cycle) + "_t"
+ cycleHead = "Cycle " + str(cycle) + " t = "
time = 0.0
for step in range(self.numberSteps + 1):
timeStampInt = int(time/self.timeUnits)
- timeStampString = repr(timeStampInt).rjust(self.timeStampWidth, '0')
+ timeStampString = repr(timeStampInt) + " years"
head += "," + cycleHead + timeStampString
time += self.deltaT
f.write('%s\n' % head)
- y = 0.0
- z = 0.0
for point in range(self.numberPoints):
- f.write(' %.12g, %.12g, %.12g' % (self.points[point], y, z))
+ f.write(' %.12g' % (self.points[point]))
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