[cig-commits] r15203 - short/3D/PyLith/trunk/tests/2d/quad4
brad at geodynamics.org
brad at geodynamics.org
Thu Jun 11 18:11:27 PDT 2009
Author: brad
Date: 2009-06-11 18:11:27 -0700 (Thu, 11 Jun 2009)
New Revision: 15203
Added:
short/3D/PyLith/trunk/tests/2d/quad4/dislocation2_soln.py
Modified:
short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation.py
short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation2.py
short/3D/PyLith/trunk/tests/2d/quad4/dislocation2.cfg
Log:
More work on 2d/quad4 full-scale tests.
Modified: short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation.py 2009-06-12 00:24:04 UTC (rev 15202)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation.py 2009-06-12 01:11:27 UTC (rev 15203)
@@ -12,7 +12,7 @@
## @file tests/2d/quad4/TestDislocation.py
##
-## @brief Test suite for testing pylith with 2-D axial extension.
+## @brief Test suite for testing pylith with fault slip.
import numpy
from TestQuad4 import TestQuad4
Modified: short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation2.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation2.py 2009-06-12 00:24:04 UTC (rev 15202)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation2.py 2009-06-12 01:11:27 UTC (rev 15203)
@@ -10,18 +10,22 @@
# ----------------------------------------------------------------------
#
-## @file tests/2d/quad4/TestDislocation2.py
+## @file tests/2d/quad4/TestDislocation.py
##
-## @brief Test suite for testing pylith with two parallel shear
-## dislocations for 2-D box.
+## @brief Test suite for testing pylith with slip on two faults.
-import unittest
+import numpy
+from TestQuad4 import TestQuad4
+from dislocation_soln import AnalyticalSoln
+from pylith.utils.VTKDataReader import has_vtk
+from pylith.utils.VTKDataReader import VTKDataReader
+from pylith.tests.Fault import check_vertex_fields
-# Local application
-from pylith.PyLithApp import PyLithApp
-class DislocationApp(PyLithApp):
+# Local version of PyLithApp
+from pylith.apps.PyLithApp import PyLithApp
+class AxialApp(PyLithApp):
def __init__(self):
- PyLithApp.__init__(self, name="dislocation2")
+ PyLithApp.__init__(self, name="dislocation")
return
@@ -31,30 +35,159 @@
Run pylith.
"""
if not "done" in dir(run_pylith):
- app = DislocationApp()
+ # Run PyLith
+ app = AxialApp()
app.run()
run_pylith.done = True
return
-class TestDislocation2(unittest.TestCase):
+class TestDislocation(TestQuad4):
"""
- Test suite for testing pylith with shear dislocation for 2-D box.
+ Test suite for testing pylith with 2-D axial extension.
"""
def setUp(self):
"""
Setup for test.
"""
+ TestQuad4.setUp(self)
+ self.mesh['nvertices'] = 81+18
+ self.nverticesO = 81
+ self.fault1Mesh = {'nvertices': 9,
+ 'spaceDim': 3,
+ 'ncells': 8,
+ 'ncorners': 2}
+ self.fault2Mesh = {'nvertices': 9,
+ 'spaceDim': 3,
+ 'ncells': 8,
+ 'ncorners': 2}
+
run_pylith()
+ self.outputRoot = "dislocation2"
+ if has_vtk():
+ self.reader = VTKDataReader()
+ self.soln = AnalyticalSoln()
+ else:
+ self.reader = None
+
return
- def test_disp(self):
+ def test_fault1_info(self):
"""
- Check displacement field.
+ Check fault information.
"""
+ if self.reader is None:
+ return
+
+ filename = "%s-fault1_info.vtk" % self.outputRoot
+ fields = ["normal_dir", "final_slip", "slip_time"]
+ check_vertex_fields(self, filename, self.fault1Mesh, fields)
+
return
+ def test_fault1_data(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault1_t0000000.vtk" % self.outputRoot
+ fields = ["cumulative_slip", "traction_change"]
+ check_vertex_fields(self, filename, self.fault1Mesh, fields)
+
+ return
+
+
+ def test_fault2_info(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault2_info.vtk" % self.outputRoot
+ fields = ["normal_dir", "final_slip", "slip_time"]
+ check_vertex_fields(self, filename, self.fault2Mesh, fields)
+
+ return
+
+
+ def test_fault2_data(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault2_t0000000.vtk" % self.outputRoot
+ fields = ["cumulative_slip", "traction_change"]
+ check_vertex_fields(self, filename, self.fault2Mesh, fields)
+
+ return
+
+
+ def calcDisplacements(self, vertices):
+ """
+ Calculate displacement field given coordinates of vertices.
+ """
+ return self.soln.displacement(vertices, self.nverticesO)
+
+
+ def calcStateVar(self, name, vertices, cells):
+ """
+ Calculate state variable.
+ """
+ ncells = self.mesh['ncells']
+ pts = numpy.zeros( (ncells, 3), dtype=numpy.float64)
+ if name == "total_strain":
+ stateVar = self.soln.strain(pts)
+ elif name == "stress":
+ stateVar = self.soln.stress(pts)
+ else:
+ raise ValueError("Unknown state variable '%s'." % name)
+
+ return stateVar
+
+
+ def calcFaultField(self, name, vertices):
+ """
+ Calculate fault info.
+ """
+
+ normalDir = (-1.0, 0.0)
+ finalSlip = -2.0
+ slipTime = 0.0
+
+ nvertices = self.faultMesh['nvertices']
+
+ if name == "normal_dir":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = normalDir[0]
+ field[:,1] = normalDir[1]
+
+ elif name == "final_slip":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = finalSlip
+
+ elif name == "slip_time":
+ field = slipTime*numpy.zeros( (nvertices, 1), dtype=numpy.float64)
+
+ elif name == "cumulative_slip":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = finalSlip
+
+ elif name == "traction_change":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = 0.0
+
+ else:
+ raise ValueError("Unknown fault field '%s'." % name)
+
+ return field
+
+
# End of file
Modified: short/3D/PyLith/trunk/tests/2d/quad4/dislocation2.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/dislocation2.cfg 2009-06-12 00:24:04 UTC (rev 15202)
+++ short/3D/PyLith/trunk/tests/2d/quad4/dislocation2.cfg 2009-06-12 01:11:27 UTC (rev 15203)
@@ -1,10 +1,13 @@
# -*- Python -*-
-[dislocation2]
+[pylithapp]
+[pylithapp.launcher] # WARNING: THIS IS NOT PORTABLE
+command = mpirun -np ${nodes}
+
# ----------------------------------------------------------------------
# journal
# ----------------------------------------------------------------------
-[dislocation2.journal.info]
+[pylithapp.journal.info]
#timedependent = 1
#implicit = 1
#petsc = 1
@@ -12,98 +15,111 @@
#meshiocubit = 1
#implicitelasticity = 1
#quadrature2d = 1
-#fiatsimplex = 1
+#fiatlagrange = 1
# ----------------------------------------------------------------------
# mesh_generator
# ----------------------------------------------------------------------
-[dislocation2.mesh_generator]
+[pylithapp.mesh_generator]
#debug = 1
-importer = pylith.meshio.MeshIOCubit
+reader = pylith.meshio.MeshIOCubit
-[dislocation2.mesh_generator.importer]
-filename = box_quad4_100m.exo
+[pylithapp.mesh_generator.reader]
+filename = mesh.exo
coordsys.space_dim = 2
# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
-[dislocation2.timedependent]
+[pylithapp.timedependent]
dimension = 2
-bc = [y_neg,y_pos]
-interfaces = [fault_pos,fault_neg]
+bc = [x_neg,x_pos]
+interfaces = [fault1,fault2]
-[dislocation2.timedependent.formulation.time_step]
+[pylithapp.timedependent.formulation.time_step]
total_time = 0.0*s
# ----------------------------------------------------------------------
# materials
# ----------------------------------------------------------------------
-[dislocation2.timedependent]
+[pylithapp.timedependent]
materials = [elastic]
materials.elastic = pylith.materials.ElasticPlaneStrain
-[dislocation2.timedependent.materials.elastic]
+[pylithapp.timedependent.materials.elastic]
label = Elastic material
id = 1
-db.iohandler.filename = matprops.spatialdb
-quadrature = pylith.feassemble.quadrature.Quadrature2D
+db_properties.iohandler.filename = matprops.spatialdb
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 2
-quadrature.cell.quad_order = 2
# ----------------------------------------------------------------------
# boundary conditions
# ----------------------------------------------------------------------
-[dislocation2.timedependent.bc.y_neg]
-fixed_dof = [0, 1]
-label = 23
-db.label = Dirichlet BC -y edge
+[pylithapp.timedependent.bc.x_pos]
+fixed_dof = [0,1]
+label = 20
+db_initial.label = Dirichlet BC +x edge
-[dislocation2.timedependent.bc.y_pos]
-fixed_dof = [0, 1]
-label = 22
-db.label = Dirichlet BC +y edge
+[pylithapp.timedependent.bc.x_neg]
+fixed_dof = [0]
+label = 21
+db_initial.label = Dirichlet BC -x edge
# ----------------------------------------------------------------------
# faults
# ----------------------------------------------------------------------
-[dislocation2.timedependent.interfaces]
+[pylithapp.timedependent.interfaces]
+fault1 = pylith.faults.FaultCohesiveKin
+fault2 = pylith.faults.FaultCohesiveKin
-[dislocation2.timedependent.interfaces.fault_pos]
+[pylithapp.timedependent.interfaces.fault1]
id = 2
-label = 11
-quadrature = pylith.feassemble.quadrature.Quadrature1Din2D
+label = 10
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 1
-quadrature.cell.quad_order = 2
-mat_db.iohandler.filename = matprops.spatialdb
-[dislocation2.timedependent.interfaces.fault_pos.eq_srcs.rupture.slip_function]
-slip.iohandler.filename = dislocation_slipLL.spatialdb
-slip_time.iohandler.filename = dislocation_sliptime.spatialdb
+[pylithapp.timedependent.interfaces.fault1.eq_srcs.rupture.slip_function]
+slip = spatialdata.spatialdb.UniformDB
+slip.label = Final slip
+slip.values = [left-lateral-slip,fault-opening]
+slip.data = [+2.0*m,0.0*m]
+slip_time = spatialdata.spatialdb.UniformDB
+slip_time.label = Slip start time
+slip_time.values = [slip-time]
+slip_time.data = [0.0*s]
-[dislocation2.timedependent.interfaces.fault_neg]
+
+[pylithapp.timedependent.interfaces.fault2]
id = 3
-label = 12
-quadrature = pylith.feassemble.quadrature.Quadrature1Din2D
+label = 11
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 1
-quadrature.cell.quad_order = 2
-mat_db.iohandler.filename = matprops.spatialdb
-[dislocation2.timedependent.interfaces.fault_neg.eq_srcs.rupture.slip_function]
-slip.iohandler.filename = dislocation_slipRL.spatialdb
-slip_time.iohandler.filename = dislocation_sliptime.spatialdb
+[pylithapp.timedependent.interfaces.fault2.eq_srcs.rupture.slip_function]
+slip = spatialdata.spatialdb.UniformDB
+slip.label = Final slip
+slip.values = [left-lateral-slip,fault-opening]
+slip.data = [-2.0*m,0.0*m]
+slip_time = spatialdata.spatialdb.UniformDB
+slip_time.label = Slip start time
+slip_time.values = [slip-time]
+slip_time.data = [0.0*s]
+
# ----------------------------------------------------------------------
# PETSc
# ----------------------------------------------------------------------
-[dislocation2.petsc]
+[pylithapp.petsc]
ksp_rtol = 1.0e-8
pc_type = asm
+# Change the preconditioner settings (must turn off
+# shift_positive_definite and turn on shift_nonzero).
+sub_pc_factor_shift_positive_definite = 0
+sub_pc_factor_shift_nonzero =
+
ksp_max_it = 100
ksp_gmres_restart = 50
#ksp_monitor = true
@@ -111,20 +127,19 @@
#log_summary = true
# start_in_debugger = true
+
# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
-# Give basename for VTK domain output of solution over domain.
-[dislocation2.problem.formulation.output.output.writer]
+[pylithapp.problem.formulation.output.output.writer]
filename = dislocation2.vtk
-# Give basename for VTK output of state variables.
-[dislocation2.timedependent.materials.elastic.output]
-cell_filter = pylith.meshio.CellFilterAvg
-writer.filename = dislocation2-statevars-elastic.vtk
+[pylithapp.timedependent.materials.elastic.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = dislocation2-elastic.vtk
-[dislocation2.timedependent.interfaces.fault_pos.output.writer]
-filename = dislocation2-fault-pos.vtk
+[pylithapp.timedependent.interfaces.fault1.output.writer]
+filename = dislocation2-fault1.vtk
-[dislocation2.timedependent.interfaces.fault_neg.output.writer]
-filename = dislocation2-fault-neg.vtk
+[pylithapp.timedependent.interfaces.fault2.output.writer]
+filename = dislocation2-fault2.vtk
Added: short/3D/PyLith/trunk/tests/2d/quad4/dislocation2_soln.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/dislocation2_soln.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/dislocation2_soln.py 2009-06-12 01:11:27 UTC (rev 15203)
@@ -0,0 +1,93 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/dislocation2_soln.py
+##
+## @brief Analytical solution to dislocation2 (two faults, rigid motion).
+
+import numpy
+
+# Physical properties
+p_density = 2500.0
+p_vs = 3000.0
+p_vp = 5291.502622129181
+
+p_mu = p_density*p_vs**2
+p_lambda = p_density*p_vp**2 - 2*p_mu
+
+# Uniform stress field (plane strain)
+sxx = 0.0
+sxy = 0.0
+syy = 0.0
+szz = p_lambda/(2*p_lambda+2*p_mu)*(sxx+syy)
+
+# Uniform strain field
+exx = 1.0/(2*p_mu) * (sxx - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+eyy = 1.0/(2*p_mu) * (syy - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+ezz = 1.0/(2*p_mu) * (szz - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+
+exy = 1.0/(2*p_mu) * (sxy)
+
+#print exx,eyy,exy,ezz
+#print -exx*p_lambda/(p_lambda+2*p_mu)
+
+# ----------------------------------------------------------------------
+class AnalyticalSoln(object):
+ """
+ Analytical solution to axial/shear displacement problem.
+ """
+
+ def __init__(self):
+ return
+
+
+ def displacement(self, locs, nlocsO):
+ """
+ Compute displacement field at locations.
+ """
+ (nlocs, dim) = locs.shape
+
+ disp = numpy.zeros( (nlocs, 3), dtype=numpy.float64)
+ maskP = locs[:,0] >= 0.0
+ maskP[nlocsO:nlocs] = False
+ maskN = numpy.bitwise_and(locs[:,0] <= 0.0, ~maskP)
+ disp[:,1] = \
+ maskN*(+1.0) + \
+ maskP*(-1.0)
+ return disp
+
+
+ def strain(self, locs):
+ """
+ Compute strain field at locations.
+ """
+ (nlocs, dim) = locs.shape
+ strain = numpy.zeros( (nlocs, 3), dtype=numpy.float64)
+ strain[:,0] = exx
+ strain[:,1] = eyy
+ strain[:,2] = exy
+ return strain
+
+
+ def stress(self, locs):
+ """
+ Compute stress field at locations.
+ """
+ (nlocs, dim) = locs.shape
+ stress = numpy.zeros( (nlocs, 3), dtype=numpy.float64)
+ stress[:,0] = sxx
+ stress[:,1] = syy
+ stress[:,2] = sxy
+ return stress
+
+
+# End of file
More information about the CIG-COMMITS
mailing list