[cig-commits] r18303 - in short/3D/PyLith/branches/v1.5-stable: libsrc/faults libsrc/friction tests/2d/frictionslide
brad at geodynamics.org
brad at geodynamics.org
Fri Apr 29 16:38:46 PDT 2011
Author: brad
Date: 2011-04-29 16:38:46 -0700 (Fri, 29 Apr 2011)
New Revision: 18303
Added:
short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/ratestate_weak.cfg
Modified:
short/3D/PyLith/branches/v1.5-stable/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.5-stable/libsrc/friction/RateStateAgeing.cc
short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/pylithapp.cfg
short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/quad4.mesh
short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/ratestate_stable.cfg
Log:
Fixed bug in rate and state friction. Limit slip rate to small infinitesimal value to keep friction from being zero.
Modified: short/3D/PyLith/branches/v1.5-stable/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/libsrc/faults/FaultCohesiveDyn.cc 2011-04-29 21:47:10 UTC (rev 18302)
+++ short/3D/PyLith/branches/v1.5-stable/libsrc/faults/FaultCohesiveDyn.cc 2011-04-29 23:38:46 UTC (rev 18303)
@@ -1922,7 +1922,7 @@
} // _sensitivitySolveLumped3D
// ----------------------------------------------------------------------
-// Constrain solution space with lumped Jacobian in 1-D.
+// Constrain solution space in 1-D.
void
pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D(double_array* dLagrangeTpdt,
const double_array& slip,
@@ -1945,7 +1945,7 @@
} // _constrainSolnSpace1D
// ----------------------------------------------------------------------
-// Constrain solution space with lumped Jacobian in 2-D.
+// Constrain solution space in 2-D.
void
pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D(double_array* dLagrangeTpdt,
const double_array& slip,
@@ -1970,7 +1970,7 @@
// Update slip based on value required to stick versus friction
const double dlp = -(tractionShearMag - frictionStress) * area *
- tractionTpdt[0] / tractionShearMag;
+ tractionTpdt[0] / tractionShearMag;
(*dLagrangeTpdt)[0] = dlp;
(*dLagrangeTpdt)[1] = 0.0;
} else {
@@ -1987,7 +1987,7 @@
} // _constrainSolnSpace2D
// ----------------------------------------------------------------------
-// Constrain solution space with lumped Jacobian in 3-D.
+// Constrain solution space in 3-D.
void
pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D(double_array* dLagrangeTpdt,
const double_array& slip,
Modified: short/3D/PyLith/branches/v1.5-stable/libsrc/friction/RateStateAgeing.cc
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/libsrc/friction/RateStateAgeing.cc 2011-04-29 21:47:10 UTC (rev 18302)
+++ short/3D/PyLith/branches/v1.5-stable/libsrc/friction/RateStateAgeing.cc 2011-04-29 23:38:46 UTC (rev 18303)
@@ -294,9 +294,14 @@
double mu_f = 0.0;
if (normalTraction <= 0.0) {
// if fault is in compression
- // Using Regularized Rate and State equation
+
+ // regularized rate and state equation
const double f0 = properties[p_coef];
+ // Since regulatized friction -> 0 as slipRate -> 0, limit slip
+ // rate to some minimum value
+ const double slipRateEff = std::max(1.0e-12, slipRate);
+
const double slipRate0 = properties[p_slipRate0];
const double a = properties[p_a];
@@ -305,7 +310,7 @@
const double b = properties[p_b];
const double bLnTerm = b * log(slipRate0 * theta / L);
const double expTerm = exp((f0 + bLnTerm)/a);
- const double sinhArg = 0.5 * slipRate / slipRate0 * expTerm;
+ const double sinhArg = 0.5 * slipRateEff / slipRate0 * expTerm;
mu_f = a * asinh(sinhArg);
friction = -mu_f * normalTraction + properties[p_cohesion];
Modified: short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/pylithapp.cfg 2011-04-29 21:47:10 UTC (rev 18302)
+++ short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/pylithapp.cfg 2011-04-29 23:38:46 UTC (rev 18303)
@@ -42,7 +42,7 @@
normalizer.length_scale = 1.0*m
normalizer.relaxation_time = 1.0*s
-bc = [top,bottom]
+bc = [xpos,xneg,ypos,yneg]
interfaces = [fault]
materials = [elastic]
@@ -50,8 +50,8 @@
solver = pylith.problems.SolverNonlinear
[pylithapp.timedependent.implicit.time_step]
-total_time = 12.0*s
-dt = 0.05*s
+total_time = 10.0*s
+dt = 0.2*s
# ----------------------------------------------------------------------
# materials
@@ -67,7 +67,6 @@
db_properties.values = [density, vp, vs]
db_properties.data = [2500.0*km/m**3, 5.1962*km/s, 3.0*km/s]
-
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 2
quadrature.cell.quad_order = 2
@@ -75,24 +74,57 @@
# ----------------------------------------------------------------------
# boundary conditions
# ----------------------------------------------------------------------
-[pylithapp.timedependent.bc.top]
+[pylithapp.timedependent.bc]
+xpos = pylith.bc.Neumann
+xneg = pylith.bc.Neumann
+ypos = pylith.bc.Neumann
+yneg = pylith.bc.DirichletBC
+
+# Dirichlet BC on -y
+[pylithapp.timedependent.bc.yneg]
bc_dof = [0, 1]
-label = top
+label = yneg
-db_rate = spatialdata.spatialdb.UniformDB
-db_rate.label = Dirichlet rate BC on +y
-db_rate.values = [displacement-rate-x,displacement-rate-y,rate-start-time]
-db_rate.data = [1.0e-6*m/s,0.0*m/s,0.0*s]
-[pylithapp.timedependent.bc.bottom]
-bc_dof = [0, 1]
-label = bottom
+# Neumann BC on +y
+[pylithapp.timedependent.bc.ypos]
+label = ypos
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+quadrature.cell.quad_order = 2
-db_rate = spatialdata.spatialdb.UniformDB
-db_rate.label = Dirichlet rate BC on -y
-db_rate.values = [displacement-rate-x,displacement-rate-y,rate-start-time]
-db_rate.data = [-1.0e-6*m/s,0.0*m/s,0.0*s]
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Neumann BC on +y
+db_initial.values = [traction-shear, traction-normal]
+db_initial.data = [-6.0*MPa, 0.0*MPa]
+
+# Neumann BC on +x
+[pylithapp.timedependent.bc.xpos]
+label = xpos
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+quadrature.cell.quad_order = 2
+
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Neumann BC on +x
+db_initial.values = [traction-shear, traction-normal]
+db_initial.data = [6.0*MPa, 0.0*MPa]
+
+
+# Neumann BC on -x
+[pylithapp.timedependent.bc.xneg]
+label = xneg
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+quadrature.cell.quad_order = 2
+
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Neumann BC on -x
+db_initial.values = [traction-shear, traction-normal]
+db_initial.data = [6.0*MPa, 0.0*MPa]
+
+
# ----------------------------------------------------------------------
# faults
# ----------------------------------------------------------------------
@@ -102,17 +134,16 @@
[pylithapp.timedependent.interfaces.fault]
id = 100
label = fault
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+quadrature.cell.quad_order = 2
db_initial_tractions = spatialdata.spatialdb.UniformDB
db_initial_tractions.label = Initial fault tractions
db_initial_tractions.values = [traction-shear,traction-normal]
-db_initial_tractions.data = [-6.1*MPa, -10.0*MPa]
+db_initial_tractions.data = [0.0*MPa, -10.0*MPa]
-quadrature.cell = pylith.feassemble.FIATLagrange
-quadrature.cell.dimension = 1
-quadrature.cell.quad_order = 2
-
# ----------------------------------------------------------------------
# PETSc
# ----------------------------------------------------------------------
@@ -123,8 +154,8 @@
sub_pc_factor_shift_type = nonzero
# KSP
-ksp_rtol = 1.0e-8
-ksp_atol = 1.0e-12
+ksp_rtol = 1.0e-12
+ksp_atol = 1.0e-15
ksp_max_it = 500
ksp_gmres_restart = 100
@@ -160,10 +191,11 @@
writer.time_format = %05.2f
[pylithapp.timedependent.interfaces.fault.output]
-vertex_data_fields=[slip,slip_rate,traction]
+vertex_data_fields=[slip,slip_rate,traction,state_variable]
writer.time_format = %05.2f
[pylithapp.timedependent.materials.elastic.output]
cell_info_fields = []
cell_data_fields = []
cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.time_format = %05.2f
Modified: short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/quad4.mesh
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/quad4.mesh 2011-04-29 21:47:10 UTC (rev 18302)
+++ short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/quad4.mesh 2011-04-29 23:38:46 UTC (rev 18303)
@@ -67,7 +67,7 @@
}
}
group = {
- name = top
+ name = ypos
type = vertices
count = 5
indices = {
@@ -75,11 +75,27 @@
}
}
group = {
- name = bottom
+ name = yneg
type = vertices
count = 5
indices = {
0 3 6 9 12
}
}
+ group = {
+ name = xpos
+ type = vertices
+ count = 3
+ indices = {
+ 12 13 14
+ }
+ }
+ group = {
+ name = xneg
+ type = vertices
+ count = 3
+ indices = {
+ 0 1 2
+ }
+ }
}
Modified: short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/ratestate_stable.cfg
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/ratestate_stable.cfg 2011-04-29 21:47:10 UTC (rev 18302)
+++ short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/ratestate_stable.cfg 2011-04-29 23:38:46 UTC (rev 18303)
@@ -29,7 +29,7 @@
friction.db_initial_state = spatialdata.spatialdb.UniformDB
friction.db_initial_state.label = Rate State Ageing State
friction.db_initial_state.values = [state-variable]
-friction.db_initial_state.data = [92.7*s]
+friction.db_initial_state.data = [1.0e+4*s]
# ----------------------------------------------------------------------
# output
Added: short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/ratestate_weak.cfg
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/ratestate_weak.cfg (rev 0)
+++ short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide/ratestate_weak.cfg 2011-04-29 23:38:46 UTC (rev 18303)
@@ -0,0 +1,45 @@
+# -*- Python -*-
+[pylithapp]
+
+# ----------------------------------------------------------------------
+# PROBLEM DESCRIPTION
+# ----------------------------------------------------------------------
+
+# Rate-state friction with weakening.
+
+# ----------------------------------------------------------------------
+# faults
+# ----------------------------------------------------------------------
+[pylithapp.timedependent.interfaces.fault]
+friction = pylith.friction.RateStateAgeing
+
+# Set the friction model parameters.
+# reference coefficient of friction: 0.6
+# reference slip rate: 1.0e-06 m/s
+# slip weakening parameter: 0.01 m
+# a: 0.008
+# b: 0.012
+# cohesion: 0 Pa
+friction.db_properties = spatialdata.spatialdb.UniformDB
+friction.db_properties.label = Rate State Ageing
+friction.db_properties.values = [reference-friction-coefficient,reference-slip-rate,characteristic-slip-distance,constitutive-parameter-a,constitutive-parameter-b,cohesion]
+friction.db_properties.data = [0.6,1.0e-6*m/s,0.01*m,0.008,0.012,0.0*Pa]
+
+# Set spatial database for the initial value of the state variable.
+friction.db_initial_state = spatialdata.spatialdb.UniformDB
+friction.db_initial_state.label = Rate State Ageing State
+friction.db_initial_state.values = [state-variable]
+friction.db_initial_state.data = [1.0e+4*s]
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+# Set filenames for output.
+[pylithapp.problem.formulation.output.output]
+writer.filename = output/ratestate_weak.vtk
+
+[pylithapp.timedependent.interfaces.fault.output]
+writer.filename = output/ratestate_weak-fault.vtk
+
+[pylithapp.timedependent.materials.elastic.output]
+writer.filename = output/ratestate_weak-statevars.vtk
More information about the CIG-COMMITS
mailing list