[cig-commits] [commit] master: Small fixes for clamped edges in cohesive cells for FaultCohesiveDyn. Added full-scale test. (11bb03b)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Sun May 25 16:02:11 PDT 2014
Repository : https://github.com/geodynamics/pylith
On branch : master
Link : https://github.com/geodynamics/pylith/compare/8727dd27d05b767d3d41da6209e9ddd606141881...11bb03bd3112722f7fb4b8f4c8beb01c5e719337
>---------------------------------------------------------------
commit 11bb03bd3112722f7fb4b8f4c8beb01c5e719337
Author: Brad Aagaard <baagaard at usgs.gov>
Date: Sun May 25 16:01:27 2014 -0700
Small fixes for clamped edges in cohesive cells for FaultCohesiveDyn. Added full-scale test.
>---------------------------------------------------------------
11bb03bd3112722f7fb4b8f4c8beb01c5e719337
libsrc/pylith/faults/FaultCohesiveDyn.cc | 8 +-
libsrc/pylith/faults/FaultCohesiveLagrange.cc | 10 +-
libsrc/pylith/friction/FrictionModel.cc | 9 ++
tests_auto/2d/tri3/Makefile.am | 2 +
...ShearDispNoSlip.py => TestShearDispFriction.py} | 43 ++++--
.../sheardispfriction.cfg} | 171 +++++++++++----------
6 files changed, 142 insertions(+), 101 deletions(-)
diff --git a/libsrc/pylith/faults/FaultCohesiveDyn.cc b/libsrc/pylith/faults/FaultCohesiveDyn.cc
index 10513c4..95be35a 100644
--- a/libsrc/pylith/faults/FaultCohesiveDyn.cc
+++ b/libsrc/pylith/faults/FaultCohesiveDyn.cc
@@ -862,6 +862,11 @@ pylith::faults::FaultCohesiveDyn::constrainSolnSpace(topology::SolutionFields* c
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
+ // Skip clamped vertices
+ if (e_lagrange < 0) {
+ continue;
+ } // if
+
// Get change in Lagrange multiplier computed from friction criterion.
const PetscInt soff = dLagrangeVisitor.sectionOffset(v_fault);
assert(spaceDim == dLagrangeVisitor.sectionDof(v_fault));
@@ -1744,8 +1749,9 @@ pylith::faults::FaultCohesiveDyn::_sensitivityUpdateJacobian(const bool negative
PetscInt* closure = NULL;
PetscInt closureSize, q = 0;
err = DMPlexGetTransitiveClosure(dmMesh, cellsCohesive[c], PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
+
// Filter out non-vertices
- for(PetscInt p = 0; p < closureSize*2; p += 2) {
+ for (PetscInt p = 0; p < closureSize*2; p += 2) {
if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
closure[q] = closure[p];
++q;
diff --git a/libsrc/pylith/faults/FaultCohesiveLagrange.cc b/libsrc/pylith/faults/FaultCohesiveLagrange.cc
index 1bd2449..a022724 100644
--- a/libsrc/pylith/faults/FaultCohesiveLagrange.cc
+++ b/libsrc/pylith/faults/FaultCohesiveLagrange.cc
@@ -1898,12 +1898,16 @@ pylith::faults::FaultCohesiveLagrange::_getJacobianSubmatrixNP(PetscMat* jacobia
int numIndicesNP = 0;
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int e_lagrange = _cohesiveVertices[iVertex].lagrange;
+ if (e_lagrange < 0) { // Ignore clamped edges.
+ continue;
+ } // if
// Compute contribution only if Lagrange constraint is local.
PetscInt goff = 0;
err = PetscSectionGetOffset(solutionGlobalSection, e_lagrange, &goff);PYLITH_CHECK_ERROR(err);
- if (goff < 0)
+ if (goff < 0) {
continue;
+ } // if
numIndicesNP += 2;
} // for
@@ -1914,6 +1918,10 @@ pylith::faults::FaultCohesiveLagrange::_getJacobianSubmatrixNP(PetscMat* jacobia
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
+ if (e_lagrange < 0) { // Ignore clamped edges.
+ continue;
+ } // if
+
// Compute contribution only if Lagrange constraint is local.
PetscInt gloff = 0;
err = PetscSectionGetOffset(solutionGlobalSection, e_lagrange, &gloff);PYLITH_CHECK_ERROR(err);
diff --git a/libsrc/pylith/friction/FrictionModel.cc b/libsrc/pylith/friction/FrictionModel.cc
index b00dd61..d84e565 100644
--- a/libsrc/pylith/friction/FrictionModel.cc
+++ b/libsrc/pylith/friction/FrictionModel.cc
@@ -28,6 +28,7 @@
#include "pylith/topology/VisitorMesh.hh" // USES VisitorMesh
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
#include "pylith/utils/array.hh" // USES scalar_array, std::vector
+#include "pylith/faults/FaultCohesiveLagrange.hh" // USES isClampedVertex()
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -194,7 +195,14 @@ pylith::friction::FrictionModel::initialize(const topology::Mesh& faultMesh,
_dbInitialState->queryVals(_metadata.dbStateVars(),
_metadata.numDBStateVars());
+ PetscDMLabel clamped = NULL;
+ PetscErrorCode err = DMPlexGetLabel(faultDMMesh, "clamped", &clamped);PYLITH_CHECK_ERROR(err);
+
for(PetscInt v = vStart; v < vEnd; ++v) {
+ if (faults::FaultCohesiveLagrange::isClampedVertex(clamped, v)) {
+ continue;
+ } // if
+
const PetscInt coff = coordsVisitor.sectionOffset(v);
assert(spaceDim == coordsVisitor.sectionDof(v));
for (PetscInt d = 0; d < spaceDim; ++d) {
@@ -225,6 +233,7 @@ pylith::friction::FrictionModel::initialize(const topology::Mesh& faultMesh,
PetscScalar* stateVarArray = stateVarVisitor.localArray();
const PetscInt off = stateVarVisitor.sectionOffset(v);
const PetscInt dof = stateVarVisitor.sectionDof(v);
+ std::cout << "v: " << v << ", dof: " << dof << ", stateVarsVetex: " << stateVarsVertex.size() << std::endl;
assert(stateVarsVertex.size() == dof);
for(PetscInt d = 0; d < dof; ++d, ++iOff) {
stateVarArray[off+d] += stateVarsVertex[iOff];
diff --git a/tests_auto/2d/tri3/Makefile.am b/tests_auto/2d/tri3/Makefile.am
index abe2d24..1f36d2e 100644
--- a/tests_auto/2d/tri3/Makefile.am
+++ b/tests_auto/2d/tri3/Makefile.am
@@ -28,6 +28,7 @@ dist_noinst_PYTHON = \
TestShearDisp.py \
TestShearDispNoSlip.py \
TestShearDispNoSlipRefine.py \
+ TestShearDispFriction.py \
sheardisp_soln.py \
sheardisp_gendb.py \
TestSlipOneFault.py \
@@ -45,6 +46,7 @@ dist_noinst_DATA = \
sheardisp.cfg \
sheardispnoslip.cfg \
sheardispnosliprefine.cfg \
+ sheardispfriction.cfg \
sliponefault.cfg \
sliptwofaults.cfg
diff --git a/tests_auto/2d/tri3/TestShearDispNoSlip.py b/tests_auto/2d/tri3/TestShearDispFriction.py
similarity index 78%
copy from tests_auto/2d/tri3/TestShearDispNoSlip.py
copy to tests_auto/2d/tri3/TestShearDispFriction.py
index 77d499f..96b08cd 100644
--- a/tests_auto/2d/tri3/TestShearDispNoSlip.py
+++ b/tests_auto/2d/tri3/TestShearDispFriction.py
@@ -16,7 +16,7 @@
# ----------------------------------------------------------------------
#
-## @file tests/2d/tri3/TestShearDispNoSlip.py
+## @file tests/2d/tri3/TestShearDispFriction.py
##
## @brief Test suite for testing pylith with 2-D shear motion with no
## fault slip.
@@ -29,7 +29,7 @@ from sheardisp_soln import AnalyticalSoln
from pylith.apps.PyLithApp import PyLithApp
class ShearApp(PyLithApp):
def __init__(self):
- PyLithApp.__init__(self, name="sheardispnoslip")
+ PyLithApp.__init__(self, name="sheardispfriction")
return
@@ -51,7 +51,7 @@ def run_pylith():
return
-class TestShearDispNoSlip(TestTri3):
+class TestShearDispFriction(TestTri3):
"""
Test suite for testing pylith with 2-D shear extension.
"""
@@ -68,7 +68,7 @@ class TestShearDispNoSlip(TestTri3):
'ncells': 2,
'ncorners': 2}
run_pylith()
- self.outputRoot = "sheardispnoslip"
+ self.outputRoot = "sheardispfriction"
self.soln = AnalyticalSoln()
return
@@ -82,7 +82,7 @@ class TestShearDispNoSlip(TestTri3):
return
filename = "%s-fault_info.h5" % self.outputRoot
- fields = ["normal_dir", "final_slip", "slip_time"]
+ fields = ["strike_dir", "normal_dir", "traction_initial","friction_coefficient","cohesion"]
from pylith.tests.Fault import check_vertex_fields
check_vertex_fields(self, filename, self.faultMesh, fields)
@@ -134,31 +134,42 @@ class TestShearDispNoSlip(TestTri3):
Calculate fault info.
"""
+ strikeDir = (1.0, 0.0)
normalDir = (0.0, 1.0)
- finalSlip = 0.0
- slipTime = 0.0
+ initialTraction = (0.0, -1.0e+7)
+ frictionCoefficient = 0.6
nvertices = self.faultMesh['nvertices']
- if name == "normal_dir":
+ if name == "strike_dir":
+ field = numpy.zeros( (1, nvertices, 2), dtype=numpy.float64)
+ field[0,:,0] = strikeDir[0]
+ field[0,:,1] = strikeDir[1]
+
+ elif name == "normal_dir":
field = numpy.zeros( (1, nvertices, 2), dtype=numpy.float64)
field[0,:,0] = normalDir[0]
field[0,:,1] = normalDir[1]
- elif name == "final_slip":
+ elif name == "traction_initial":
field = numpy.zeros( (1, nvertices, 2), dtype=numpy.float64)
- field[0,:,0] = finalSlip
+ field[0,:,0] = initialTraction[0]
+ field[0,:,1] = initialTraction[1]
+
+ elif name == "friction_coefficient":
+ field = numpy.zeros( (1, nvertices, 1), dtype=numpy.float64)
+ field[0,:,0] = frictionCoefficient
- elif name == "slip_time":
- field = slipTime*numpy.zeros( (1, nvertices, 1), dtype=numpy.float64)
+ elif name == "cohesion":
+ field = numpy.zeros( (1, nvertices, 1), dtype=numpy.float64)
elif name == "slip":
field = numpy.zeros( (1, nvertices, 2), dtype=numpy.float64)
- field[0,:,0] = finalSlip
- elif name == "traction_change":
+ elif name == "traction":
field = numpy.zeros( (1, nvertices, 2), dtype=numpy.float64)
- field[0,:,0] = 0.0
+ field[0,:,0] = 1.0e+6
+ field[0,:,1] = -1.0e+7
else:
raise ValueError("Unknown fault field '%s'." % name)
@@ -175,7 +186,7 @@ class TestShearDispNoSlip(TestTri3):
# ----------------------------------------------------------------------
if __name__ == '__main__':
import unittest
- from TestShearDispNoSlip import TestShearDispNoSlip as Tester
+ from TestShearDispFriction import TestShearDispFriction as Tester
suite = unittest.TestSuite()
suite.addTest(unittest.makeSuite(Tester))
diff --git a/tests_auto/2d/quad4/friction_opening.cfg b/tests_auto/2d/tri3/sheardispfriction.cfg
similarity index 54%
copy from tests_auto/2d/quad4/friction_opening.cfg
copy to tests_auto/2d/tri3/sheardispfriction.cfg
index 02cdc5e..c5e5775 100644
--- a/tests_auto/2d/quad4/friction_opening.cfg
+++ b/tests_auto/2d/tri3/sheardispfriction.cfg
@@ -1,105 +1,108 @@
-[friction_opening]
+[sheardispfriction]
+nodes = 1
-[friction_opening.launcher] # WARNING: THIS IS NOT PORTABLE
+[sheardispfriction.launcher] # WARNING: THIS IS NOT PORTABLE
command = mpirun -np ${nodes}
# ----------------------------------------------------------------------
# journal
# ----------------------------------------------------------------------
-# The settings below turn on journal info for the specified components.
-# If you want less output to stdout, you can turn these off.
-[friction_opening.journal.info]
-#timedependent = 1
-#implicit = 1
-#petsc = 1
-#solvernonlinear = 1
-#meshioascii = 1
-#homogeneous = 1
-#elasticityimplicit = 1
-#fiatlagrange = 1
-#quadrature1d = 1
-#faultcohesivedyn = 1
+[sheardispfriction.journal.info]
+sheardispfriction = 1
+timedependent = 1
+implicit = 1
+petsc = 1
+solvernonlinear = 1
+meshimporter = 1
+meshiocubit = 1
+implicitelasticity = 1
+quadrature2d = 1
+faultcohesivekin = 1
+fiatsimplex = 1
# ----------------------------------------------------------------------
# mesh_generator
# ----------------------------------------------------------------------
-# The settings below control the mesh generation (importing mesh info).
-# Turn on debugging output for mesh generation.
-[friction_opening.mesh_generator]
-#debug = 1
+[sheardispfriction.mesh_generator]
reader = pylith.meshio.MeshIOCubit
+reorder_mesh = True
-# This component specification means we are using PyLith ASCII format,
-# and we then specify the filename and number of space dimensions for
-# the mesh.
-[friction_opening.mesh_generator.reader]
+[sheardispfriction.mesh_generator.reader]
filename = mesh.exo
-use_nodeset_names = False
coordsys.space_dim = 2
# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
-[friction_opening.timedependent]
+[sheardispfriction.timedependent]
dimension = 2
-formulation = pylith.problems.Implicit
+normalizer.length_scale = 1.0*km
formulation.solver = pylith.problems.SolverNonlinear
-bc = [x_neg,x_pos]
-bc.x_pos=pylith.bc.DirichletBoundary
-
-interfaces = [fault]
-
-
-[friction_opening.timedependent.formulation.time_step]
+[sheardispfriction.timedependent.formulation.time_step]
total_time = 0.0*s
-dt = 10.0*year
-
# ----------------------------------------------------------------------
# materials
# ----------------------------------------------------------------------
-[friction_opening.timedependent.materials]
-material = pylith.materials.ElasticPlaneStrain
+[sheardispfriction.timedependent]
+materials = [elastic]
+materials.elastic = pylith.materials.ElasticPlaneStrain
-[friction_opening.timedependent.materials.material]
+[sheardispfriction.timedependent.materials.elastic]
label = Elastic material
id = 1
-
db_properties.label = Elastic properties
db_properties.iohandler.filename = matprops.spatialdb
-
-quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2
# ----------------------------------------------------------------------
# boundary conditions
# ----------------------------------------------------------------------
-[friction_opening.timedependent.bc.x_neg]
-bc_dof = [0,1]
-label = 21
-
-[friction_opening.timedependent.bc.x_pos]
-bc_dof = [0,1]
-label = 20
+[sheardispfriction.timedependent]
+bc = [x_neg,x_pos,y_neg,y_pos]
-db_initial = spatialdata.spatialdb.UniformDB
+[sheardispfriction.timedependent.bc.x_pos]
+bc_dof = [1]
+label = edge_xpos
+db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Dirichlet BC +x edge
-db_initial.values = [displacement-x,displacement-y]
-db_initial.data = [1.0*m,0.0*m]
+db_initial.iohandler.filename = shear_disp.spatialdb
+
+[sheardispfriction.timedependent.bc.x_neg]
+bc_dof = [1]
+label = edge_xneg
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC -x edge
+db_initial.iohandler.filename = shear_disp.spatialdb
+
+[sheardispfriction.timedependent.bc.y_pos]
+bc_dof = [0]
+label = edge_ypos
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC +y edge
+db_initial.iohandler.filename = shear_disp.spatialdb
+
+[sheardispfriction.timedependent.bc.y_neg]
+bc_dof = [0]
+label = edge_yneg
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC -y edge
+db_initial.iohandler.filename = shear_disp.spatialdb
# ----------------------------------------------------------------------
# faults
# ----------------------------------------------------------------------
-# Provide information on the fault (interface).
-[friction_opening.timedependent.interfaces]
-
-fault = pylith.faults.FaultCohesiveDyn
-
-[friction_opening.timedependent.interfaces.fault]
-label = 10
+[sheardispfriction.timedependent]
+interfaces = [fault]
+interfaces.fault = pylith.faults.FaultCohesiveDyn
-quadrature.cell = pylith.feassemble.FIATLagrange
+[sheardispfriction.timedependent.interfaces.fault]
+id = 100
+label = fault_y
+edge = fault_y_edge
+quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 1
friction.label = Static friction
@@ -110,39 +113,37 @@ friction.db_properties.data = [0.6,0.0*Pa]
traction_perturbation = pylith.faults.TractPerturbation
-[friction_opening.timedependent.interfaces.fault.traction_perturbation]
+[sheardispfriction.timedependent.interfaces.fault.traction_perturbation]
db_initial = spatialdata.spatialdb.UniformDB
db_initial.label = Initial fault tractions
db_initial.values = [traction-shear,traction-normal]
-db_initial.data = [0.0*Pa, -1.0*MPa]
+db_initial.data = [0.0*Pa, -10.0*MPa]
# ----------------------------------------------------------------------
# PETSc
# ----------------------------------------------------------------------
-[friction_opening.petsc]
+[sheardispfriction.petsc]
pc_type = asm
# Change the preconditioner settings.
-sub_pc_factor_shift_type = nonzero
+sub_pc_factor_shift_type = none
-ksp_rtol = 1.0e-8
-ksp_atol = 1.0e-12
+ksp_rtol = 1.0e-20
+ksp_atol = 1.0e-10
ksp_max_it = 100
ksp_gmres_restart = 50
-#ksp_monitor = true
-#ksp_view = true
-#ksp_converged_reason = true
-snes_rtol = 1.0e-8
-snes_atol = 1.0e-9
+snes_rtol = 1.0e-20
+snes_atol = 1.0e-8
snes_max_it = 100
-#snes_monitor = true
-#snes_ls_monitor = true
-#snes_view = true
-#snes_converged_reason = true
-#log_summary = true
-#start_in_debugger = true
+ksp_monitor = true
+#ksp_view = true
+ksp_converged_reason = true
+
+snes_monitor = true
+#snes_view = true
+snes_converged_reason = true
# Friction sensitivity solve used to compute the increment in slip
# associated with changes in the Lagrange multiplier imposed by the
@@ -155,19 +156,23 @@ friction_ksp_gmres_restart = 30
#friction_ksp_monitor = true
#friction_ksp_view = true
#friction_ksp_converged_reason = true
+
+#log_summary = true
+#start_in_debugger = true
+
# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
-[friction_opening.problem.formulation.output.output]
+[sheardispfriction.problem.formulation.output.output]
writer = pylith.meshio.DataWriterHDF5
-writer.filename = friction_opening.h5
+writer.filename = sheardispfriction.h5
-[friction_opening.timedependent.interfaces.fault.output]
+[sheardispfriction.timedependent.interfaces.fault.output]
writer = pylith.meshio.DataWriterHDF5
-writer.filename = friction_opening-fault.h5
-vertex_info_fields = [strike_dir,normal_dir,traction_initial_value,friction_coefficient,cohesion]
+writer.filename = sheardispfriction-fault.h5
-[friction_opening.timedependent.materials.material.output]
+[sheardispfriction.timedependent.materials.elastic.output]
cell_filter = pylith.meshio.CellFilterAvg
writer = pylith.meshio.DataWriterHDF5
-writer.filename = friction_opening-elastic.h5
+writer.filename = sheardispfriction-elastic.h5
+
More information about the CIG-COMMITS
mailing list