[cig-commits] r19973 - in short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott: . figures
willic3 at geodynamics.org
willic3 at geodynamics.org
Mon Apr 23 16:55:37 PDT 2012
Author: willic3
Date: 2012-04-23 16:55:37 -0700 (Mon, 23 Apr 2012)
New Revision: 19973
Modified:
short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/calc_analytic.cfg
short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/calc_analytic.py
short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/figures/plot_displ_profiles.py
Log:
Fixed some problems in calc_analytic.py.
Changed number of time steps per cycle and explicitly put in filenames.
Started creating script to plot displacement profiles but didn't finish.
-This line, and those below, will be ignored--
M savageprescott/figures/plot_displ_profiles.py
M savageprescott/calc_analytic.py
M savageprescott/calc_analytic.cfg
Modified: short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/calc_analytic.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/calc_analytic.cfg 2012-04-23 23:13:21 UTC (rev 19972)
+++ short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/calc_analytic.cfg 2012-04-23 23:55:37 UTC (rev 19973)
@@ -10,8 +10,13 @@
plate_velocity = 2.0*cm/year
number_cycles = 10
-number_steps = 40
+number_steps = 20
number_terms = 100
number_points = 400
delta_x = 5.0*km
x_epsilon = 0.01*m
+
+output_disp = True
+disp_filename = analytical/analytic_disp.txt
+output_vel = True
+vel_filename = analytical/analytic_vel.txt
Modified: short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/calc_analytic.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/calc_analytic.py 2012-04-23 23:13:21 UTC (rev 19972)
+++ short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/calc_analytic.py 2012-04-23 23:55:37 UTC (rev 19973)
@@ -101,13 +101,13 @@
outputVel = pyre.inventory.bool("output_vel", default=True)
outputVel.meta['tip'] = "Output velocity file?"
- dispfilename = pyre.inventory.str("disp_filename",
- default="output/analytic_disp.txt")
- dispfilename.meta['tip'] = "Filename for displacement output."
+ dispFilename = pyre.inventory.str("disp_filename",
+ default="analytical/analytic_disp.txt")
+ dispFilename.meta['tip'] = "Filename for displacement output."
- velfilename = pyre.inventory.str("vel_filename",
- default="output/analytic_vel.txt")
- velfilename.meta['tip'] = "Filename for velocity output."
+ velFilename = pyre.inventory.str("vel_filename",
+ default="analytical/analytic_vel.txt")
+ velFilename.meta['tip'] = "Filename for velocity output."
# PUBLIC METHODS /////////////////////////////////////////////////////
@@ -148,11 +148,11 @@
self.outputDisp = self.inventory.outputDisp
self.outputVel = self.inventory.outputVel
- self.dispfilename = self.inventory.dispfilename
- self.velfilename = self.inventory.velfilename
+ self.dispFilename = self.inventory.dispFilename
+ self.velFilename = self.inventory.velFilename
- self.deltaT = self.recurrenceTime/self.numberSteps
- self.tauFac = 0.5*self.shearModulus/self.viscosity
+ self.deltaT = self.recurrenceTime / self.numberSteps
+ self.tauFac = 0.5 * self.shearModulus / self.viscosity
self.tau0 = self.recurrenceTime * self.tauFac
return
@@ -163,6 +163,8 @@
Create array of points for output along with series terms
for each point.
"""
+ print "Generating points and point coefficients:"
+
self.points = numpy.zeros(self.numberPoints, dtype=numpy.float64)
self.pointCoeff = numpy.zeros((self.numberPoints, self.numberTerms),
dtype=numpy.float64)
@@ -188,6 +190,8 @@
"""
Compute transient solution.
"""
+ print "Generating solution:"
+
solutionU2 = numpy.zeros((self.numberCycles,
self.numberSteps + 1,
self.numberPoints),
@@ -212,6 +216,7 @@
for cycle in range(self.numberCycles):
+ print " Working on cycle # %d:" % cycle
time = cycle * self.numberSteps * deltaT
tau = time * tauFac
if cycle > 0:
@@ -293,7 +298,7 @@
Computes viscoelastic solution for times greater than the recurrence time.
"""
tau0 = self.tau0
- velocity = self.velocity
+ velocity = self.velocity.value
recurrenceTime = self.recurrenceTime.value
solutionU = numpy.zeros(self.numberPoints, dtype=numpy.float64)
@@ -334,10 +339,12 @@
"""
if solutionType == "displacement":
- filename = self.dispfilename
+ print "Writing displacements:"
+ filename = self.dispFilename
solution = self.solutionUTot
else:
- filename = self.velfilename
+ print "Writing velocities:"
+ filename = self.velFilename
solution = self.solutionVTot
from pyre.units.time import year
@@ -346,10 +353,10 @@
f.write("dt = %f*year\n" % (self.deltaT.value/year.value,))
f.write("dx = %s\n" % self.deltaX)
f.write("ncycles = %d\n" % self.numberCycles)
- f.write("nsteps = %d\n" % (self.numberSteps+1,))
+ f.write("nsteps = %d\n" % (self.numberSteps + 1,))
f.write("npoints = %d\n" % self.numberPoints)
data = solution.reshape( (self.numberCycles,
- (self.numberSteps+1)*self.numberPoints) )
+ (self.numberSteps + 1) * self.numberPoints) )
numpy.savetxt(f, data, fmt="%14.6e")
f.close()
Modified: short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/figures/plot_displ_profiles.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/figures/plot_displ_profiles.py 2012-04-23 23:13:21 UTC (rev 19972)
+++ short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott/figures/plot_displ_profiles.py 2012-04-23 23:55:37 UTC (rev 19973)
@@ -10,49 +10,91 @@
# ======================================================================
#
-cycles = [2,9]
+analyticalFile = "../analytical/savpres_displ.csv"
+numericalFiles = ["../output/sp_hex8/sp_hex8-groundsurf.h5",
+ "../output/sp_tet4/sp_tet4-groundsurf.h5"]
+# The indices for analytical solution need changing for time step size of 10.
+anlRefCol = 370
+anlOutputCols = [372, 380, 390, 400, 408]
+# End change.
+numRefStep = 180
+numOutputSteps = [181, 185, 190, 195, 199]
+numSteps = 5
+numPlots = 2
+anlHeaderLines = 5
+
style = "lightbg"
-fileSuffix = ".eps"
plotSize = "poster"
-numericalFileRoot = "../results/spbm_hex8_graded_20km/spbm_hex8_graded_c"
-analyticalFileRoot = "../utils/savpres_displ_c"
-outputFileRoot = "spbm_hex8_graded_c"
-coseismicDispl = 400.0
-elasticThick = 40.0
+outputFiles = ["spbm_hex8.eps", "spbm_tet4.eps"]
+coseismicDispl = 4.0
+elasticThick = 4.0e4
# ======================================================================
import numpy
import pylab
+import tables
from Figure import Figure
# ----------------------------------------------------------------------
class ProfileSet(object):
"""
- Set of profiles.
+ Get profiles for analytical solution.
"""
- def __init__(self, filename):
- self.data = numpy.loadtxt(filename, comments="#")
- ncols = self.data.shape[1]
- self.data[:,1:ncols+1] /= coseismicDispl # Coseismic displacement
- self.data[:,0] /= elasticThick # Thickness of elastic layer
+ def __init__(self, filename, solnType):
+ if (solnType == "analytical"):
+ # Analytical solution stuff is completely broken for the new analytical
+ # output. I was just starting to fix it.
+ f = open(filename, 'r')
+ for line in range(anlHeaderLines):
+
+ soln = numpy.loadtxt(filename, comments="#", dtype=numpy.float64)
+ self.xCoord = soln[:,0]/elasticThick
+ numPoints = self.xCoord.shape[0]
+ refDispl = soln[:,anlRefCol]
+ self.data = (soln[:,anlOutputCols] - refDispl)/coseismicDispl
+ else:
+ # Open file and get coordinates
+ h5 = tables.openFile(filename, "r")
+ coords = h5.root.geometry.vertices[:]
+
+ # Get coordinate indices for y == 0.0 and x >= 0.0
+ yProf = coords[:,1] == 0.0
+ posVals = coords[:,0] >= 0.0
+ posProf = numpy.logical_and(yProf, posVals)
+ posProfInds = numpy.nonzero(posProf)
+
+ # Get solution, and then remove profile indices that give negative displ.
+ soln = h5.root.vertex_fields.displacement[:]
+ solnTest = soln[numRefStep, posProfInds, 1]
+ negInd = numpy.nonzero(solnTest < 0.0)[1][0]
+ profileInds = numpy.delete(posProfInds, negInd)
+
+ # Get test profile coordinates and then sort key.
+ profCoords = coords[profileInds, 0]
+ coordSort = numpy.argsort(profCoords)
+ profInds = profileInds[coordSort]
+
+ # Get coordinates and solution.
+ self.xCoord = coords[profInds, 0]/elasticThick
+ numPoints = self.xCoord.shape[0]
+ refDispl = soln[numRefStep, profInds, 1]
+ refDisplReshape = numpy.column_stack((refDispl, refDispl, refDispl,
+ refDispl, refDispl))
+ solnSteps = soln[numOutputSteps, :, 1]
+ self.data = (solnSteps[:, profInds].transpose() - refDisplReshape)/coseismicDispl
+
return
def plot(self, linestyle, linewidth):
data = self.data
- ncols = data.shape[1]
- if ncols == 42:
- indicesTime = [2, 10, 20, 30, 38]
- elif ncols == 21:
- indicesTime = [1, 5, 10, 15, 19]
+ xCoord = self.xCoord
colorOrder = ('orange', 'blue', 'red', 'green', 'purple')
- i = 0
- for index in indicesTime:
- h = pylab.plot(data[:,0], data[:,1+index], linestyle, linewidth=linewidth,
- color=colorOrder[i])
- i += 1
+ for step in range(numSteps):
+ h = pylab.plot(xCoord, data[:,step], linestyle, linewidth=linewidth,
+ color=colorOrder[step])
return h
@@ -62,7 +104,7 @@
Figure with displacement time histories for a site.
"""
- def __init__(self, cycle):
+ def __init__(self, plotNum):
if plotSize == "poster":
fontsize = 14
elif plotSize == "presentation":
@@ -72,7 +114,8 @@
else:
raise ValueError("Unknown plotSize '%s'." % plotSize)
Figure.__init__(self, color=style, fontsize=fontsize)
- self.cycle = cycle
+ self.plotNum = plotNum
+
return
@@ -87,12 +130,11 @@
# Plot profiles
self.axes(1, 1, 1, 1)
- analyticFile = analyticalFileRoot + repr(self.cycle) + ".txt"
- analytic = ProfileSet(analyticFile)
+ analytic = ProfileSet(analyticalFile, "analytical")
ah = analytic.plot(linestyle="-", linewidth=1)
pylab.hold(True)
- simulationFile = numericalFileRoot + repr(self.cycle) + ".txt"
- simulation = ProfileSet(simulationFile)
+ simulationFile = numericalFiles[self.plotNum]
+ simulation = ProfileSet(simulationFile, "numerical")
sh = simulation.plot(linestyle="--", linewidth=1)
pylab.hold(False)
@@ -110,7 +152,7 @@
def save(self):
pylab.figure(self.handle.number)
- outputFile = outputFileRoot + repr(self.cycle) + fileSuffix
+ outputFile = outputFiles[self.plotNum]
pylab.savefig(outputFile)
return
@@ -133,12 +175,12 @@
# ----------------------------------------------------------------------
if __name__ == "__main__":
- # import pdb
- # pdb.set_trace()
+ import pdb
+ pdb.set_trace()
figures = []
- for cycle in cycles:
- figure = Profiles(cycle)
+ for plotNum in range(numPlots):
+ figure = Profiles(plotNum)
figure.plot()
figures.append(figure)
More information about the CIG-COMMITS
mailing list