[cig-commits] r19304 - in short/3D/PyLith/branches/v1.6-stable: examples/3d/hex8 libsrc/pylith/faults libsrc/pylith/friction tests/2d/frictionslide
brad at geodynamics.org
brad at geodynamics.org
Mon Dec 19 17:44:08 PST 2011
Author: brad
Date: 2011-12-19 17:44:07 -0800 (Mon, 19 Dec 2011)
New Revision: 19304
Added:
short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/test.cfg
Modified:
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc
short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py
short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg
short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate.cfg
short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_stable.cfg
short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_weak.cfg
Log:
More work on rate state friction.
Added: short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/test.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/test.cfg (rev 0)
+++ short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/test.cfg 2011-12-20 01:44:07 UTC (rev 19304)
@@ -0,0 +1,215 @@
+# -*- Python -*-
+[pylithapp]
+
+# ----------------------------------------------------------------------
+# PROBLEM DESCRIPTION
+# ----------------------------------------------------------------------
+
+#
+# This problem demonstrates the use of rate-and-state friction for a
+# quasi-static problem, using the aging law for evolution of the state
+# variable.
+#
+# The problem is similar to example 13, except that a different
+# friction model is used. The axial deformation and time step have
+# also been changed. The model is run for 200 years.
+
+# ----------------------------------------------------------------------
+# RUNNING THE SIMULATON
+# ----------------------------------------------------------------------
+
+# This is not a self-contained simulation configuration file. This
+# file only specifies parameters specific to tutorial step14.
+# The general parameters are specificed in the pylithapp.cfg
+# file which PyLith reads by default.
+#
+# To run the simulation:
+# pylith step14.cfg
+#
+# Output will be directed to the directory output.
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[pylithapp.timedependent]
+# Set bc to an array of 3 boundary conditions: 'x_pos', 'x_neg', and 'z_neg'.
+bc = [x_pos,x_neg,z_neg]
+
+# Set interfaces to an array of 1 fault: 'fault'.
+interfaces = [fault]
+
+[pylithapp.timedependent.implicit]
+# Set the output to an array of 2 output managers.
+# We will output the solution over the domain and the ground surface.
+output = [domain,subdomain]
+
+# Set subdomain component to OutputSolnSubset (boundary of the domain).
+output.subdomain = pylith.meshio.OutputSolnSubset
+
+# Fault friction is a nonlinear problem so we need to use the nonlinear
+# solver.
+solver = pylith.problems.SolverNonlinear
+
+# Change the total simulation time to 200 years, and use a constant time
+# step size of 5 years.
+[pylithapp.timedependent.implicit.time_step]
+total_time = 200.0*year
+dt = 0.5*year
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+# Set the parameters for boundary conditions applied on the
+# +x, -x, and -z faces of the box.
+#
+# On the -x and +x faces, we fix the x degrees of freedom and apply
+# velocities in the y-direction. We fix the z degree of freedom on the
+# bottom (-z) face.
+#
+# We use a UniformDB to apply the displacements and velocities, while
+# retaining the default ZeroDispDB for zero x-displacements on the -x
+# face and and zero z-displacements on -z.
+#
+# Note that since the fault cuts through the base of the model (z_neg),
+# we can only constrain the portion of the bottom boundary that does not
+# include the fault. A nodeset named 'face_zneg_nofault' has been defined
+# in Cubit for this purpose.
+#
+
+# The label corresponds to the name of the nodeset in CUBIT.
+
+# +x face -- Dirichlet
+[pylithapp.timedependent.bc.x_pos]
+bc_dof = [0,1]
+label = face_xpos
+
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Dirichlet BC on +x
+db_initial.values = [displacement-x,displacement-y]
+db_initial.data = [-2.0*m,0.0*m]
+
+db_rate = spatialdata.spatialdb.UniformDB
+db_rate.label = Dirichlet rate BC on +x
+db_rate.values = [displacement-rate-x,displacement-rate-y,rate-start-time]
+db_rate.data = [0.0*cm/year,1.0*cm/year,0.0*year]
+
+# -x face
+[pylithapp.timedependent.bc.x_neg]
+bc_dof = [0, 1]
+label = face_xneg
+db_initial.label = Dirichlet BC on -x
+db_rate = spatialdata.spatialdb.UniformDB
+db_rate.label = Dirichlet rate BC on -x
+db_rate.values = [displacement-rate-x,displacement-rate-y,rate-start-time]
+db_rate.data = [0.0*cm/year,-1.0*cm/year,0.0*year]
+
+# -z face
+[pylithapp.timedependent.bc.z_neg]
+bc_dof = [2]
+label = face_zneg_nofault
+db_initial.label = Dirichlet BC on -z
+
+# ----------------------------------------------------------------------
+# faults
+# ----------------------------------------------------------------------
+[pylithapp.timedependent.interfaces]
+# Change fault to dynamic fault interface.
+fault = pylith.faults.FaultCohesiveDyn
+
+[pylithapp.timedependent.interfaces.fault]
+# The label corresponds to the name of the nodeset in CUBIT.
+label = fault
+zero_tolerance = 1.0e-15
+
+# Use the rate-and-state aging friction model.
+friction = pylith.friction.RateStateAgeing
+friction.label = Rate and state
+friction.min_slip_rate = 1.0e-10
+
+# We must define the quadrature information for fault cells.
+# The fault cells are 2D (surface).
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 2
+
+# Set rate-and-state parameters using a uniform DB. Set the parameters as
+# follows:
+# reference coefficient of friction: 0.4
+# reference slip rate: 1.0e-9 m/s
+# characteristic slip distance: 0.05 m
+# a: 0.01
+# b: 0.08
+# cohesion: 0 Pa
+friction.db_properties = spatialdata.spatialdb.UniformDB
+friction.db_properties.label = Rate Stete 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.4, 1.0e-9*m/s, 0.05*m, 0.01, 0.08, 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]
+# theta_ss = characteristic_slip_dist / reference_slip_rate
+friction.db_initial_state.data = [5.0e+07*s]
+
+# ----------------------------------------------------------------------
+# PETSc settings
+# ----------------------------------------------------------------------
+# NOTE: There are additional settings specific to fault friction.
+[pylithapp.petsc]
+
+# Friction sensitivity solve used to compute the increment in slip
+# associated with changes in the Lagrange multiplier imposed by the
+# fault constitutive model.
+friction_pc_type = asm
+friction_sub_pc_factor_shift_type = nonzero
+friction_ksp_max_it = 25
+friction_ksp_gmres_restart = 30
+# Uncomment to view details of friction sensitivity solve.
+#friction_ksp_monitor = true
+#friction_ksp_view = true
+friction_ksp_converged_reason = true
+
+# Reduce convergence tolerances.
+ksp_rtol = 1.0e-18
+ksp_atol = 1.0e-18
+
+snes_rtol = 1.0e-16
+snes_atol = 1.0e-16
+
+snes_max_it = 200
+
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[pylithapp.problem.formulation.output.domain]
+output_freq = time_step
+time_step = 10.0*year
+writer = pylith.meshio.DataWriterHDF5Mesh
+writer.filename = output/test.h5
+
+[pylithapp.problem.formulation.output.subdomain]
+label = face_zpos
+writer = pylith.meshio.DataWriterHDF5SubMesh
+writer.filename = output/test-groundsurf.h5
+
+[pylithapp.problem.interfaces.fault.output]
+writer = pylith.meshio.DataWriterHDF5SubSubMesh
+writer.filename = output/test-fault.h5
+vertex_data_fields = [slip, slip_rate, traction, state_variable]
+
+[pylithapp.timedependent.materials.upper_crust.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+output_freq = time_step
+time_step = 20.0*year
+writer = pylith.meshio.DataWriterHDF5Mesh
+writer.filename = output/test-upper_crust.h5
+
+# Give basename for VTK output of lower_crust state variables.
+[pylithapp.timedependent.materials.lower_crust.output]
+# Average values over quadrature points.
+cell_filter = pylith.meshio.CellFilterAvgMesh
+output_freq = time_step
+time_step = 20.0*year
+writer = pylith.meshio.DataWriterHDF5Mesh
+writer.filename = output/test-lower_crust.h5
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-12-19 22:22:19 UTC (rev 19303)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-12-20 01:44:07 UTC (rev 19304)
@@ -353,9 +353,11 @@
residualSection->getFiberDimension(v_positive));
residualSection->updateAddPoint(v_positive, &residualVertexP[0]);
+#if 0 // TEST
assert(residualVertexL.size() ==
residualSection->getFiberDimension(v_lagrange));
residualSection->updateAddPoint(v_lagrange, &residualVertexL[0]);
+#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
@@ -526,7 +528,7 @@
double_array dDispRelVertex(spaceDim);
// Get sections
- double_array slipVertex(spaceDim);
+ double_array slipTpdtVertex(spaceDim);
const ALE::Obj<RealSection>& dispRelSection =
_fields->get("relative disp").section();
assert(!dispRelSection.isNull());
@@ -629,12 +631,12 @@
// Compute slip, slip rate, and Lagrange multiplier at time t+dt
// in fault coordinate system.
- slipVertex = 0.0;
+ slipTpdtVertex = 0.0;
slipRateVertex = 0.0;
tractionTpdtVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
- slipVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ slipTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
(dispTVertexP[jDim] + dispIncrVertexP[jDim]
- dispTVertexN[jDim] - dispIncrVertexN[jDim]);
slipRateVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
@@ -646,43 +648,43 @@
slipRateVertex[iDim] = 0.0;
} // if
} // for
- if (fabs(slipVertex[indexN]) < _zeroTolerance) {
- slipVertex[indexN] = 0.0;
+ if (fabs(slipTpdtVertex[indexN]) < _zeroTolerance) {
+ slipTpdtVertex[indexN] = 0.0;
} // if
double dSlipVertexNormal = 0.0;
double dTractionTpdtVertexNormal = 0.0;
- if (slipVertex[indexN]*tractionTpdtVertex[indexN] < 0.0) {
+ if (slipTpdtVertex[indexN]*tractionTpdtVertex[indexN] < 0.0) {
#if 0 // DEBUGGING
std::cout << "STEP 1 CORRECTING NONPHYSICAL SLIP/TRACTIONS"
<< ", v_fault: " << v_fault
- << ", slipNormal: " << slipVertex[indexN]
+ << ", slipNormal: " << slipTpdtVertex[indexN]
<< ", tractionNormal: " << tractionTpdtVertex[indexN]
<< std::endl;
#endif
// Don't know what behavior is appropriate so set smaller of
// traction and slip to zero (should be appropriate if problem
// is nondimensionalized correctly).
- if (fabs(slipVertex[indexN]) > fabs(tractionTpdtVertex[indexN])) {
+ if (fabs(slipTpdtVertex[indexN]) > fabs(tractionTpdtVertex[indexN])) {
// slip is bigger, so force normal traction back to zero
dTractionTpdtVertexNormal = -tractionTpdtVertex[indexN];
tractionTpdtVertex[indexN] = 0.0;
} else {
// traction is bigger, so force slip back to zero
- dSlipVertexNormal = -slipVertex[indexN];
- slipVertex[indexN] = 0.0;
+ dSlipVertexNormal = -slipTpdtVertex[indexN];
+ slipTpdtVertex[indexN] = 0.0;
} // if/else
} // if
- if (slipVertex[indexN] < 0.0) {
+ if (slipTpdtVertex[indexN] < 0.0) {
#if 0 // DEBUGGING
std::cout << "STEP 1 CORRECTING INTERPENETRATION"
<< ", v_fault: " << v_fault
- << ", slipNormal: " << slipVertex[indexN]
+ << ", slipNormal: " << slipTpdtVertex[indexN]
<< ", tractionNormal: " << tractionTpdtVertex[indexN]
<< std::endl;
#endif
- dSlipVertexNormal = -slipVertex[indexN];
- slipVertex[indexN] = 0.0;
+ dSlipVertexNormal = -slipTpdtVertex[indexN];
+ slipTpdtVertex[indexN] = 0.0;
} // if
// Step 2: Apply friction criterion to trial solution to get
@@ -698,7 +700,7 @@
const bool iterating = true; // Iterating to get friction
CALL_MEMBER_FN(*this,
constrainSolnSpaceFn)(&dLagrangeTpdtVertex,
- slipVertex, slipRateVertex,
+ slipTpdtVertex, slipRateVertex,
tractionTpdtVertex,
iterating);
@@ -720,7 +722,7 @@
std::cout << "v_fault: " << v_fault;
std::cout << ", slipVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << slipVertex[iDim];
+ std::cout << " " << slipTpdtVertex[iDim];
std::cout << ", slipRateVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << slipRateVertex[iDim];
@@ -794,7 +796,8 @@
// Lagrange multipliers) and Step 3 (slip associated with change in
// Lagrange multipliers).
- double_array dSlipVertex(spaceDim);
+ double_array slipTVertex(spaceDim);
+ double_array dSlipTpdtVertex(spaceDim);
double_array dispRelVertex(spaceDim);
const ALE::Obj<RealSection>& sensDispRelSection =
@@ -869,66 +872,92 @@
// interpenetration with tension or opening with compression).
// Compute slip, change in slip, and tractions in fault coordinates.
- dSlipVertex = 0.0;
- slipVertex = 0.0;
+ slipTVertex = 0.0;
+ slipTpdtVertex = 0.0;
+ dSlipTpdtVertex = 0.0;
tractionTpdtVertex = 0.0;
dTractionTpdtVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
- dSlipVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
- sensDispRelVertex[jDim];
- slipVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ slipTVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ (dispTVertexP[jDim] - dispTVertexN[jDim]);
+ slipTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
(dispTVertexP[jDim] - dispTVertexN[jDim] +
dispIncrVertexP[jDim] - dispIncrVertexN[jDim]);
+ dSlipTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ sensDispRelVertex[jDim];
tractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
(lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim]);
dTractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
dLagrangeTpdtVertex[jDim];
} // for
} // for
- if (fabs(slipVertex[indexN]) < _zeroTolerance) {
- slipVertex[indexN] = 0.0;
+ if (fabs(slipTpdtVertex[indexN]) < _zeroTolerance) {
+ slipTpdtVertex[indexN] = 0.0;
} // if
- if (fabs(dSlipVertex[indexN]) < _zeroTolerance) {
- dSlipVertex[indexN] = 0.0;
+ if (fabs(dSlipTpdtVertex[indexN]) < _zeroTolerance) {
+ dSlipTpdtVertex[indexN] = 0.0;
} // if
- if ((slipVertex[indexN] + dSlipVertex[indexN]) *
+ if ((slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN]) *
(tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN])
< 0.0) {
#if 0 // DEBUGGING
std::cout << "STEP 4a CORRECTING NONPHYSICAL SLIP/TRACTIONS"
<< ", v_fault: " << v_fault
- << ", slipNormal: " << slipVertex[indexN] + dSlipVertex[indexN]
+ << ", slipNormal: " << slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN]
<< ", tractionNormal: " << tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN]
<< std::endl;
#endif
// Don't know what behavior is appropriate so set smaller of
// traction and slip to zero (should be appropriate if problem
// is nondimensionalized correctly).
- if (fabs(slipVertex[indexN] + dSlipVertex[indexN]) >
+ if (fabs(slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN]) >
fabs(tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN])) {
// slip is bigger, so force normal traction back to zero
dTractionTpdtVertex[indexN] = -tractionTpdtVertex[indexN];
} else {
// traction is bigger, so force slip back to zero
- dSlipVertex[indexN] = -slipVertex[indexN];
+ dSlipTpdtVertex[indexN] = -slipTpdtVertex[indexN];
} // if/else
} // if
// Do not allow fault interpenetration.
- if (slipVertex[indexN] + dSlipVertex[indexN] < 0.0) {
+ if (slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN] < 0.0) {
#if 0 // DEBUGGING
std::cout << "STEP 4a CORRECTING INTERPENETATION"
<< ", v_fault: " << v_fault
- << ", slipNormal: " << slipVertex[indexN] + dSlipVertex[indexN]
+ << ", slipNormal: " << slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN]
<< std::endl;
#endif
- dSlipVertex[indexN] = -slipVertex[indexN];
+ dSlipTpdtVertex[indexN] = -slipTpdtVertex[indexN];
} // if
+ // Prevent over-correction in slip resulting in backslip.
+ double slipDot = 0.0;
+ double tractionSlipDot = 0.0;
+
+ // :TODO:
+ //
+ // Need slipTVertex, slipTpdtVertex, dSlipTpdtVertex
+ // slipDot = (slipTpdtVertex - slipTVertex) dot dSlipVertex
+ //
+ // if slipDot < 0, dSlipVertex = -0.5*(slipTpdtVertex - slipTVertex)
+
+ for (int iDim=0; iDim < spaceDim-1; ++iDim) { // :TODO: Optimize this
+ slipDot += (slipTpdtVertex[iDim] - slipTVertex[iDim]) * (slipTpdtVertex[iDim] + dSlipTpdtVertex[iDim] - slipTVertex[iDim]);
+ tractionSlipDot += (tractionTpdtVertex[iDim] + dTractionTpdtVertex[iDim])
+ * (slipTpdtVertex[iDim] + dSlipTpdtVertex[iDim]);
+ } // for
+ if (slipDot < 0.0 && tractionSlipDot < 0.0) {
+ dTractionTpdtVertex *= 0.5; // Use bisection as guess for traction
+ for (int iDim=0; iDim < spaceDim-1; ++iDim) {
+ dSlipTpdtVertex[iDim] *= -0.5*(slipTpdtVertex[iDim] - slipTVertex[iDim]);
+ } // for
+ } // if
+
// Update current estimate of slip from t to t+dt.
- slipVertex += dSlipVertex;
+ slipTpdtVertex += dSlipTpdtVertex;
// Compute relative displacement from slip.
dispRelVertex = 0.0;
@@ -937,9 +966,9 @@
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
dispRelVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
- slipVertex[jDim];
+ slipTpdtVertex[jDim];
dDispRelVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
- dSlipVertex[jDim];
+ dSlipTpdtVertex[jDim];
dLagrangeTpdtVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
dTractionTpdtVertex[jDim];
} // for
@@ -948,7 +977,7 @@
dDispTIncrVertexP[iDim] = +0.5*dDispRelVertex[iDim];
} // for
-#if 0 // debugging
+#if 1 // debugging
std::cout << "v_fault: " << v_fault;
std::cout << ", tractionTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
@@ -956,26 +985,29 @@
std::cout << ", dTractionTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << dTractionTpdtVertex[iDim];
- std::cout << ", dLagrangeTpdtVertex: ";
+ //std::cout << ", dLagrangeTpdtVertex: ";
+ //for (int iDim=0; iDim < spaceDim; ++iDim)
+ // std::cout << " " << dLagrangeTpdtVertex[iDim];
+ std::cout << ", slipTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dLagrangeTpdtVertex[iDim];
- std::cout << ", slipVertex: ";
+ std::cout << " " << slipTpdtVertex[iDim];
+ std::cout << ", dSlipTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << slipVertex[iDim];
- std::cout << ", dispRelVertex: ";
- for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dispRelVertex[iDim];
- std::cout << ", dDispRelVertex: ";
- for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dDispRelVertex[iDim];
+ std::cout << " " << dSlipTpdtVertex[iDim];
+ //std::cout << ", dDispRelVertex: ";
+ //for (int iDim=0; iDim < spaceDim; ++iDim)
+ // std::cout << " " << dDispRelVertex[iDim];
+ std::cout << ", slipDot: " << slipDot
+ << ", tractionSlipDot: " << tractionSlipDot;
std::cout << std::endl;
#endif
+#if 0 // TEST
// Set change in relative displacement.
assert(dispRelVertex.size() ==
dispRelSection->getFiberDimension(v_fault));
dispRelSection->updatePoint(v_fault, &dispRelVertex[0]);
-
+#endif
// Update Lagrange multiplier increment.
assert(dLagrangeTpdtVertex.size() ==
dispIncrSection->getFiberDimension(v_lagrange));
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc 2011-12-19 22:22:19 UTC (rev 19303)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc 2011-12-20 01:44:07 UTC (rev 19304)
@@ -72,7 +72,7 @@
const int numDBStateVars = 1;
const char* dbStateVars[1] = {
- "state-variable"
+ "state-variable",
};
} // _RateStateAgeing
@@ -255,9 +255,7 @@
const int numDBValues = dbValues.size();
assert(_RateStateAgeing::numDBStateVars == numDBValues);
- const double stateVariable = dbValues[db_state];
-
- stateValues[s_state] = stateVariable;
+ stateValues[s_state] = dbValues[db_state];
} // _dbToStateVars
// ----------------------------------------------------------------------
@@ -311,6 +309,7 @@
if (normalTraction <= 0.0) {
// if fault is in compression
+#if 0
// regularized rate and state equation
const double f0 = properties[p_coef];
@@ -322,14 +321,57 @@
const double a = properties[p_a];
const double theta = stateVars[s_state];
+
const double L = properties[p_L];
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 * slipRateEff / slipRate0 * expTerm;
mu_f = a * asinh(sinhArg);
friction = -mu_f * normalTraction + properties[p_cohesion];
+
+ std::cout << "slip: " << slip
+ << ", slipRate: " << slipRate
+ << ", stateVar: " << theta
+ << ", bLnTermL: " << bLnTerm
+ << ", expTerm: " << expTerm
+ << ", sinhArg: " << sinhArg
+ << ", mu_f: " << mu_f
+ << std::endl;
+#else
+
+ const double slipRateLinear = _minSlipRate;
+ const double slipRateFactor = 1.0e-3;
+
+ const double f0 = properties[p_coef];
+ const double a = properties[p_a];
+ const double b = properties[p_b];
+ const double L = properties[p_L];
+ const double slipRate0 = properties[p_slipRate0];
+ const double theta = stateVars[s_state];
+
+ if (slipRate > slipRateLinear) {
+ mu_f = f0 + a*log(slipRate / slipRate0) + b*log(slipRate0*theta/L);
+ } else if (slipRate > slipRateFactor*slipRateLinear) {
+ mu_f = f0 + a*log(slipRateLinear / slipRate0) + b*log(slipRate0*theta/L) -
+ a*(1.0-slipRateFactor) *
+ (1.0 - slipRate/slipRateLinear) / (1.0 - slipRateFactor);
+ } else {
+ mu_f = f0 + a*log(slipRateLinear / slipRate0) + b*log(slipRate0*theta/L) -
+ a*(1.0-slipRateFactor);
+ } // else
+
+ friction = -mu_f * normalTraction + properties[p_cohesion];
+
+ std::cout << "slip: " << slip
+ << ", slipRate: " << slipRate
+ << ", stateVar: " << theta
+ << ", mu_f: " << mu_f
+ << std::endl;
+
+#endif
} // if
PetscLogFlops(11);
@@ -374,7 +416,7 @@
// Since regulatized friction -> 0 as slipRate -> 0, limit slip
// rate to some minimum value
- const double slipRateEff = std::max(_minSlipRate, slipRate);
+ const double slipRateEff = std::max(_minSlipRate, slipRate);
const double dt = _dt;
const double thetaTVertex = stateVars[s_state];
Modified: short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py 2011-12-19 22:22:19 UTC (rev 19303)
+++ short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py 2011-12-20 01:44:07 UTC (rev 19304)
@@ -1,6 +1,6 @@
#!/usr/bin/env python
-sim = "ratestate_stable"
+sim = "ratestate_weak"
# ======================================================================
import tables
Modified: short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg 2011-12-19 22:22:19 UTC (rev 19303)
+++ short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg 2011-12-20 01:44:07 UTC (rev 19304)
@@ -47,7 +47,7 @@
[pylithapp.timedependent.implicit.time_step]
total_time = 14.0*s
-dt = 0.05*s
+dt = 0.1*s
# ----------------------------------------------------------------------
# materials
Modified: short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate.cfg 2011-12-19 22:22:19 UTC (rev 19303)
+++ short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate.cfg 2011-12-20 01:44:07 UTC (rev 19304)
@@ -46,7 +46,7 @@
db_initial_tractions.values = [traction-shear,traction-normal]
db_initial_tractions.data = [0.0*MPa, -10.0*MPa]
-zero_tolerance = 1.0e-14
+zero_tolerance = 1.0e-16
# ----------------------------------------------------------------------
# PETSc
Modified: short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_stable.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_stable.cfg 2011-12-19 22:22:19 UTC (rev 19303)
+++ short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_stable.cfg 2011-12-20 01:44:07 UTC (rev 19304)
@@ -13,6 +13,7 @@
[pylithapp.timedependent.interfaces.fault]
friction = pylith.friction.RateStateAgeing
friction.label = Rate state ageing
+friction.min_slip_rate = 1.0e-8
# Set the friction model parameters.
# reference coefficient of friction: 0.6
@@ -30,7 +31,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 = [2.0e+0*s]
+friction.db_initial_state.data = [1.5*s]
# ----------------------------------------------------------------------
# output
Modified: short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_weak.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_weak.cfg 2011-12-19 22:22:19 UTC (rev 19303)
+++ short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_weak.cfg 2011-12-20 01:44:07 UTC (rev 19304)
@@ -13,7 +13,7 @@
[pylithapp.timedependent.interfaces.fault]
friction = pylith.friction.RateStateAgeing
friction.label = Rate state ageing
-friction.min_slip_rate = 1.0e-14
+friction.min_slip_rate = 1.0e-8
# Set the friction model parameters.
# reference coefficient of friction: 0.6
More information about the CIG-COMMITS
mailing list