[cig-commits] r20263 - in short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns: reverse strikeslip
willic3 at geodynamics.org
willic3 at geodynamics.org
Thu May 31 22:19:15 PDT 2012
Author: willic3
Date: 2012-05-31 22:19:14 -0700 (Thu, 31 May 2012)
New Revision: 20263
Added:
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/make-plot.sh
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_results.py
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/make-plot.sh
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_results.py
Modified:
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/README
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/README
Log:
Added simple python scripts to plot inversion results, and added stuff to
README's.
Modified: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/README
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/README 2012-06-01 03:48:15 UTC (rev 20262)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/README 2012-06-01 05:19:14 UTC (rev 20263)
@@ -69,9 +69,30 @@
Step 3. Invert for coseismic slip
- ADD STUFF HERE
+ Now that we have the forward simulation (data) and Green's functions, we
+ can perform an inversion for the fault slip. We have written a very
+ simple Python script that performs a miminim moment solution using a
+ range of penalty parameters (contained in penalty_params.txt). The Python
+ code is in simple_invert.py, and a script to run it is contained in
+ run-inv.sh. You can run the inversion by simply typing:
+ ./run-inv.sh
+ This will generate a set of predicted slip values, contained in
+ strikeslip_inversion.txt.
+
+Step 4. Plot the true and predicted slip values.
+
+ To see how the predicted solution compares with the applied slip, we
+ use a simple Python script (plot_results.py) that makes use of
+ matplotlib. You can run this script by simply typing:
+ ./make-plot.sh
+
+ This will bring up a matplotlib plot window with the true solution (shown
+ in black) and the predicted solution corresponding to different values
+ of the penalty parameter.
+
+
Suggestions variations
The list below includes some suggested modifications to the problem
@@ -86,3 +107,8 @@
* Adjust the location and number of the observation points.
* Increase/decrease the resolution of the mesh.
+
+ * Try a different penalty function for the inversion. One simple option
+ would be a 1D Laplacian approximation (sets of values of [-1, 2, -1])
+ centered along the diagonal.
+
Added: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/make-plot.sh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/make-plot.sh (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/make-plot.sh 2012-06-01 05:19:14 UTC (rev 20263)
@@ -0,0 +1,2 @@
+#!/bin/bash
+./plot_results.py --solution_file=output/eqsim-fault.h5 --predicted_file=reverse_inversion.txt
Property changes on: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/make-plot.sh
___________________________________________________________________
Name: svn:executable
+ *
Added: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_results.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_results.py (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_results.py 2012-06-01 05:19:14 UTC (rev 20263)
@@ -0,0 +1,85 @@
+#!/usr/bin/env python
+"""
+This script generates a set of line plots comparing the predicted solution to
+the true solution.
+"""
+
+# The code requires the numpy, h5py, and matplotlib packages.
+import numpy
+import h5py
+import matplotlib.pyplot as pyplot
+# import pdb
+
+# Define line colors and other parameters.
+solnColor = "black"
+predColors = ["cyan", "red", "blue", "green", "orange", "purple", "yellow"]
+style = "lightbg"
+fontsize = 8
+
+
+def getSolution():
+ """
+ Function to get applied slip and fault coordinates.
+ """
+
+ # Open solution file and get slip and coordinates.
+ solution = h5py.File(solutionFile, "r")
+ solC = solution['geometry/vertices'][:]
+ solV = solution['vertex_fields/slip'][:,:,0].flatten()
+
+ # Sort by y-coordinate.
+ solInds = numpy.argsort(solC[:,1])
+ solCSort = solC[solInds,:]
+ solVSort = solV[solInds]
+ # Reverse sign, since this is right-lateral slip.
+ solVSort *= -1.0
+ solution.close()
+
+ return (solCSort, solVSort)
+
+
+# pdb.set_trace()
+
+# The main part of the code is below.
+# Get command-line arguments.
+from optparse import OptionParser
+parser = OptionParser()
+parser.add_option("-i", "--solution_file", action="store", type="string",
+ dest="solution_file",
+ help="HDF5 file with applied fault slip")
+parser.add_option("-r", "--predicted_file", action="store", type="string",
+ dest="predicted_file", help="file with predicted solutions")
+
+(options, args) = parser.parse_args()
+
+if (not options.solution_file):
+ parser.error("solution input file must be specified")
+
+if (not options.predicted_file):
+ parser.error("predicted input file must be specified")
+
+solutionFile = options.solution_file
+predictedFile = options.predicted_file
+
+# Get applied fault slip.
+(solCoords, solVals) = getSolution()
+
+# Get predicted fault slip.
+predicted = numpy.loadtxt(predictedFile, skiprows=1, dtype=numpy.float64)
+predCoords = predicted[:,0:2]
+predSlip = predicted[:, 2:]
+
+# Generate figure, starting with true solution.
+pyplot.plot(solCoords[:,1], solVals, color=solnColor)
+pyplot.xlabel("Y-distance (m)")
+pyplot.ylabel("Reverse slip (m)")
+
+# Loop over predicted results, using a different line color for each.
+numInv = predSlip.shape[1]
+
+for inversion in range(numInv):
+ pyplot.plot(predCoords[:,1], predSlip[:,inversion],
+ color=predColors[inversion])
+
+
+pyplot.show()
Property changes on: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_results.py
___________________________________________________________________
Name: svn:executable
+ *
Modified: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/README
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/README 2012-06-01 03:48:15 UTC (rev 20262)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/README 2012-06-01 05:19:14 UTC (rev 20263)
@@ -69,9 +69,30 @@
Step 3. Invert for coseismic slip
- ADD STUFF HERE
+ Now that we have the forward simulation (data) and Green's functions, we
+ can perform an inversion for the fault slip. We have written a very
+ simple Python script that performs a miminim moment solution using a
+ range of penalty parameters (contained in penalty_params.txt). The Python
+ code is in simple_invert.py, and a script to run it is contained in
+ run-inv.sh. You can run the inversion by simply typing:
+ ./run-inv.sh
+ This will generate a set of predicted slip values, contained in
+ strikeslip_inversion.txt.
+
+Step 4. Plot the true and predicted slip values.
+
+ To see how the predicted solution compares with the applied slip, we
+ use a simple Python script (plot_results.py) that makes use of
+ matplotlib. You can run this script by simply typing:
+ ./make-plot.sh
+
+ This will bring up a matplotlib plot window with the true solution (shown
+ in black) and the predicted solution corresponding to different values
+ of the penalty parameter.
+
+
Suggestions variations
The list below includes some suggested modifications to the problem
@@ -86,3 +107,7 @@
* Adjust the location and number of the observation points.
* Increase/decrease the resolution of the mesh.
+
+ * Try a different penalty function for the inversion. One simple option
+ would be a 1D Laplacian approximation (sets of values of [-1, 2, -1])
+ centered along the diagonal.
Added: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/make-plot.sh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/make-plot.sh (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/make-plot.sh 2012-06-01 05:19:14 UTC (rev 20263)
@@ -0,0 +1,2 @@
+#!/bin/bash
+./plot_results.py --solution_file=output/eqsim-fault.h5 --predicted_file=strikeslip_inversion.txt
Property changes on: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/make-plot.sh
___________________________________________________________________
Name: svn:executable
+ *
Added: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_results.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_results.py (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_results.py 2012-06-01 05:19:14 UTC (rev 20263)
@@ -0,0 +1,85 @@
+#!/usr/bin/env python
+"""
+This script generates a set of line plots comparing the predicted solution to
+the true solution.
+"""
+
+# The code requires the numpy, h5py, and matplotlib packages.
+import numpy
+import h5py
+import matplotlib.pyplot as pyplot
+# import pdb
+
+# Define line colors and other parameters.
+solnColor = "black"
+predColors = ["cyan", "red", "blue", "green", "orange", "purple", "yellow"]
+style = "lightbg"
+fontsize = 8
+
+
+def getSolution():
+ """
+ Function to get applied slip and fault coordinates.
+ """
+
+ # Open solution file and get slip and coordinates.
+ solution = h5py.File(solutionFile, "r")
+ solC = solution['geometry/vertices'][:]
+ solV = solution['vertex_fields/slip'][:,:,0].flatten()
+
+ # Sort by y-coordinate.
+ solInds = numpy.argsort(solC[:,1])
+ solCSort = solC[solInds,:]
+ solVSort = solV[solInds]
+ # Reverse sign, since this is right-lateral slip.
+ solVSort *= -1.0
+ solution.close()
+
+ return (solCSort, solVSort)
+
+
+# pdb.set_trace()
+
+# The main part of the code is below.
+# Get command-line arguments.
+from optparse import OptionParser
+parser = OptionParser()
+parser.add_option("-i", "--solution_file", action="store", type="string",
+ dest="solution_file",
+ help="HDF5 file with applied fault slip")
+parser.add_option("-r", "--predicted_file", action="store", type="string",
+ dest="predicted_file", help="file with predicted solutions")
+
+(options, args) = parser.parse_args()
+
+if (not options.solution_file):
+ parser.error("solution input file must be specified")
+
+if (not options.predicted_file):
+ parser.error("predicted input file must be specified")
+
+solutionFile = options.solution_file
+predictedFile = options.predicted_file
+
+# Get applied fault slip.
+(solCoords, solVals) = getSolution()
+
+# Get predicted fault slip.
+predicted = numpy.loadtxt(predictedFile, skiprows=1, dtype=numpy.float64)
+predCoords = predicted[:,0:2]
+predSlip = predicted[:, 2:]
+
+# Generate figure, starting with true solution.
+pyplot.plot(solCoords[:,1], solVals, color=solnColor)
+pyplot.xlabel("Y-distance (m)")
+pyplot.ylabel("Right-lateral slip (m)")
+
+# Loop over predicted results, using a different line color for each.
+numInv = predSlip.shape[1]
+
+for inversion in range(numInv):
+ pyplot.plot(predCoords[:,1], predSlip[:,inversion],
+ color=predColors[inversion])
+
+
+pyplot.show()
Property changes on: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_results.py
___________________________________________________________________
Name: svn:executable
+ *
More information about the CIG-COMMITS
mailing list