[cig-commits] r8620 - short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils

[brad at geodynamics.org]

Author: brad
Date: 2007-12-07 17:15:14 -0800 (Fri, 07 Dec 2007)
New Revision: 8620

Added:
   short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/calc_analytic.m
   short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/compare_solns.py
Modified:
   short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/vtktoh5.py
Log:
Added scripts for driving analytic solution calculation and doing comparisons.

Added: short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/calc_analytic.m
===================================================================
--- short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/calc_analytic.m	2007-12-08 00:44:46 UTC (rev 8619)
+++ short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/calc_analytic.m	2007-12-08 01:15:14 UTC (rev 8620)
@@ -0,0 +1,7 @@
+addpath('../../okada');
+scenario = 'tet4_0250m';
+
+in = sprintf('../analytic/output/%s_points.in', scenario);
+out = sprintf('../analytic/output/%s_points.disp', scenario);
+
+StrikeSlip_ng(in, out);

Added: short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/compare_solns.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/compare_solns.py	2007-12-08 00:44:46 UTC (rev 8619)
+++ short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/compare_solns.py	2007-12-08 01:15:14 UTC (rev 8620)
@@ -0,0 +1,119 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+#                           Brad T. Aagaard
+#                        U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @brief Python script to compare benchmark solutions.
+
+import tables
+
+shape = "tet4"
+res = 500
+
+# ----------------------------------------------------------------------
+def getMap(solnfile):
+    """
+    Get element interpolation mapping information.
+    """
+    vertices = solnfile.root.geometry.vertices[:]
+    (nvertices, spaceDim) = vertices.shape
+    cells = solnfile.root.topology.cells[:]
+    (ncells, ncorners) = cells.shape
+
+    if spaceDim == 3 and ncorners == 4:
+        from benchmark.discretize.Tet4 import Tet4
+        cell = Tet4()
+    else:
+        raise ValueError("Unknown finite-element in dimension %d with %d " \
+                         "corners." % (spaceDim, ncorners))
+    return cell.calcMap(vertices, cells)
+
+
+# ----------------------------------------------------------------------
+def main_load_disp(app):
+    """
+    Local replacement for CompareApp::main().
+
+    Load displacement fields.
+    """
+    
+    field = "displacements"
+    filename = "results/strikeslip_%s_%04dm.h5" % (shape, res)
+    
+    projection = app.projection
+    projection.open()
+
+    # PyLith ---------------------------------
+    app._info.log("Projecting PyLith solution...")
+    solnfile = tables.openFile(filename, 'r')
+    solnmap = getMap(solnfile)
+    for tstep in [0]:
+        projection.project(solnfile, solnmap, "pylith_1_0", tstep, field)
+    solnfile.close()
+
+    # Analytic -------------------------------
+    app._info.log("Copying analytic solution...")
+    filename = "analytic/output/%s_%04dm.h5" % (shape, res)
+    solnfile = tables.openFile(filename, 'r')
+    for tstep in [0]:
+        projection.copy_projection(solnfile, "analytic", tstep, field)
+    solnfile.close()
+
+    # ----------------------------------------
+    projection.close()
+
+    return
+
+# ----------------------------------------------------------------------
+def main_cmp_disp(app):
+    """
+    Local replacement for CompareApp::main().
+
+    Compare displacement fields.
+    """
+    
+    field = "displacements"
+    
+    projection = app.projection
+    
+    projection.open()
+
+    # Code versus analytic -------------------
+    tstep = 0
+    print "Time Step: %d" % tstep
+
+    err = projection.difference(tstep, field, "pylith_1_0", "analytic")
+    print "PyLith-1.0, error wrt analytic: %12.4e" % err
+    
+    # ----------------------------------------
+    projection.close()
+
+    return
+
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+
+    from pkg_resources import require
+    require("merlin>=1.0")
+    require("pythia>=0.8")
+
+    import journal
+    journal.info("compareapp").activate()
+
+    from benchmark.CompareApp import CompareApp
+    app = CompareApp()
+    app.mainfn = main_load_disp
+    app.run()
+    app.mainfn = main_cmp_disp
+    app.run()
+
+
+# End of file


Property changes on: short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/compare_solns.py
___________________________________________________________________
Name: svn:executable
   + *

Modified: short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/vtktoh5.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/vtktoh5.py	2007-12-08 00:44:46 UTC (rev 8619)
+++ short/3D/PyLith/benchmarks/trunk/quasistatic/strikeslipnog/utils/vtktoh5.py	2007-12-08 01:15:14 UTC (rev 8620)
@@ -41,10 +41,6 @@
 
 
 h5 = tables.openFile(filenameOut, "w")
-h5.root._v_attrs.space_dim = spaceDim
-h5.root._v_attrs.num_corners = ncorners
-h5.root._v_attrs.num_vertices = nvertices
-h5.root._v_attrs.cell_shape = shape
 
 h5.createGroup("/", "topology")
 h5.createArray("/topology", "cells", cellsVtk)