[cig-commits] r19485 - in short/3D/PyLith/trunk: . examples/3d/hex8 libsrc/pylith/faults libsrc/pylith/friction libsrc/pylith/meshio libsrc/pylith/problems libsrc/pylith/topology modulesrc/friction pylith/apps pylith/feassemble pylith/friction tests/2d/frictionslide tests/3d tests/3d/cyclicfriction unittests/libtests/faults
brad at geodynamics.org
brad at geodynamics.org
Thu Jan 26 17:09:28 PST 2012
Author: brad
Date: 2012-01-26 17:09:27 -0800 (Thu, 26 Jan 2012)
New Revision: 19485
Added:
short/3D/PyLith/trunk/examples/3d/hex8/test.cfg
short/3D/PyLith/trunk/tests/3d/refine/
Modified:
short/3D/PyLith/trunk/README
short/3D/PyLith/trunk/examples/3d/hex8/step13.cfg
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/TopologyOps.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/MeshRefiner.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/RefineEdges2.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/RefineFace4Edges2.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/RefineUniform.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/RefineVol8Face4Edges2.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/ReverseCuthillMcKee.cc
short/3D/PyLith/trunk/modulesrc/friction/SlipWeakening.i
short/3D/PyLith/trunk/pylith/apps/PyLithApp.py
short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
short/3D/PyLith/trunk/pylith/friction/RateStateAgeing.py
short/3D/PyLith/trunk/pylith/friction/SlipWeakening.py
short/3D/PyLith/trunk/pylith/friction/TimeWeakening.py
short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg
short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate.cfg
short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate_stable.cfg
short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate_weak.cfg
short/3D/PyLith/trunk/tests/2d/frictionslide/tension.cfg
short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
Log:
Merge from stable.
Modified: short/3D/PyLith/trunk/README
===================================================================
--- short/3D/PyLith/trunk/README 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/README 2012-01-27 01:09:27 UTC (rev 19485)
@@ -1,4 +1,4 @@
-We are pleased to announce release of PyLith version 1.6.2
+We are pleased to announce release of PyLith version 1.6.3
Please submit bug reports via the World Wide Web at:
http://geodynamics.org/roundup
Modified: short/3D/PyLith/trunk/examples/3d/hex8/step13.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/step13.cfg 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/examples/3d/hex8/step13.cfg 2012-01-27 01:09:27 UTC (rev 19485)
@@ -8,9 +8,9 @@
#
# This problem demonstrates the use of slip-weakening fault friction for a
# quasi-static problem.
+#
# The problem is similar to example 12, except that a different friction
-# model is used. The model is run for 200 years. Slip begins to occur
-# at about 80 years, and continues in each step after that.
+# model is used. The model is run for 200 years.
# ----------------------------------------------------------------------
# RUNNING THE SIMULATON
@@ -121,6 +121,9 @@
# Use the slip-weakening friction model.
friction = pylith.friction.SlipWeakening
friction.label = Slip weakening
+# Force healing after each time step, so weakening is confined to each
+# time step and does not carry over into subsequent time steps.
+friction.force_healing = True
# We must define the quadrature information for fault cells.
# The fault cells are 2D (surface).
Copied: short/3D/PyLith/trunk/examples/3d/hex8/test.cfg (from rev 19484, short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/test.cfg)
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/test.cfg (rev 0)
+++ short/3D/PyLith/trunk/examples/3d/hex8/test.cfg 2012-01-27 01:09:27 UTC (rev 19485)
@@ -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/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -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);
@@ -530,7 +532,7 @@
scalar_array dDispRelVertex(spaceDim);
// Get sections
- scalar_array slipVertex(spaceDim);
+ scalar_array slipTpdtVertex(spaceDim);
const ALE::Obj<RealSection>& dispRelSection =
_fields->get("relative disp").section();
assert(!dispRelSection.isNull());
@@ -633,12 +635,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] *
@@ -650,43 +652,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
PylithScalar dSlipVertexNormal = 0.0;
PylithScalar 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
@@ -702,7 +704,7 @@
const bool iterating = true; // Iterating to get friction
CALL_MEMBER_FN(*this,
constrainSolnSpaceFn)(&dLagrangeTpdtVertex,
- t, slipVertex, slipRateVertex,
+ t, slipTpdtVertex, slipRateVertex,
tractionTpdtVertex,
iterating);
@@ -721,9 +723,10 @@
} // for
#if 0 // debugging
- std::cout << "slipVertex: ";
+ 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];
@@ -797,7 +800,8 @@
// Lagrange multipliers) and Step 3 (slip associated with change in
// Lagrange multipliers).
- scalar_array dSlipVertex(spaceDim);
+ scalar_array slipTVertex(spaceDim);
+ scalar_array dSlipTpdtVertex(spaceDim);
scalar_array dispRelVertex(spaceDim);
const ALE::Obj<RealSection>& sensDispRelSection =
@@ -812,6 +816,14 @@
dLagrangeTpdtSection->restrictPoint(v_fault, &dLagrangeTpdtVertex[0],
dLagrangeTpdtVertex.size());
+ // If no change in the Lagrange multiplier computed from friction
+ // criterion, there are no updates, so continue.
+ double dLagrangeTpdtMag = 0.0;
+ for (int i=0; i < spaceDim; ++i)
+ dLagrangeTpdtMag += dLagrangeTpdtVertex[i]*dLagrangeTpdtVertex[i];
+ if (0.0 == dLagrangeTpdtMag)
+ continue;
+
// Get change in relative displacement from sensitivity solve.
assert(spaceDim == sensDispRelSection->getFiberDimension(v_fault));
const PylithScalar* sensDispRelVertex =
@@ -864,66 +876,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;
@@ -932,9 +970,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
@@ -943,7 +981,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)
@@ -951,26 +989,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));
@@ -2257,6 +2298,7 @@
// if in compression and no opening
const PylithScalar frictionStress =
_friction->calcFriction(t, slipShearMag, slipRateMag, tractionNormal);
+
if (tractionShearMag > frictionStress || (iterating && slipRateMag > 0.0)) {
// traction is limited by friction, so have sliding OR
// friction exceeds traction due to overshoot in slip
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/TopologyOps.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/TopologyOps.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/TopologyOps.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -169,7 +169,7 @@
}
// This only works for uninterpolated meshes
- assert((mesh->depth() == 1) || (mesh->depth() == -1));
+ assert((mesh->depth() == 1) || (mesh->depth() == 0));
ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> sV(std::max(1, sieve->getMaxSupportSize()));
ALE::ISieveVisitor::PointRetriever<SieveMesh::sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
for(PointSet::const_iterator fv_iter = fvBegin; fv_iter != fvEnd; ++fv_iter) {
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -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 PylithScalar stateVariable = dbValues[db_state];
-
- stateValues[s_state] = stateVariable;
+ stateValues[s_state] = dbValues[db_state];
} // _dbToStateVars
// ----------------------------------------------------------------------
@@ -312,6 +310,7 @@
if (normalTraction <= 0.0) {
// if fault is in compression
+#if 0
// regularized rate and state equation
const PylithScalar f0 = properties[p_coef];
@@ -323,14 +322,57 @@
const PylithScalar a = properties[p_a];
const PylithScalar theta = stateVars[s_state];
+
const PylithScalar L = properties[p_L];
const PylithScalar b = properties[p_b];
const PylithScalar bLnTerm = b * log(slipRate0 * theta / L);
+
const PylithScalar expTerm = exp((f0 + bLnTerm)/a);
const PylithScalar 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);
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -112,7 +112,8 @@
_SlipWeakening::stateVars,
_SlipWeakening::numStateVars,
_SlipWeakening::dbStateVars,
- _SlipWeakening::numDBStateVars))
+ _SlipWeakening::numDBStateVars)),
+ _forceHealing(false)
{ // constructor
} // constructor
@@ -125,6 +126,14 @@
// ----------------------------------------------------------------------
// Compute properties from values in spatial database.
void
+pylith::friction::SlipWeakening::forceHealing(const bool flag)
+{ // forceHealing
+ _forceHealing = flag;
+} // forceHealing
+
+// ----------------------------------------------------------------------
+// Compute properties from values in spatial database.
+void
pylith::friction::SlipWeakening::_dbToProperties(
PylithScalar* const propValues,
const scalar_array& dbValues) const
@@ -274,11 +283,14 @@
PylithScalar mu_f = 0.0;
if (normalTraction <= 0.0) {
// if fault is in compression
- if (stateVars[s_slipCum] < properties[p_d0]) {
+ const double slipPrev = stateVars[s_slipPrev];
+ const double slipCum = stateVars[s_slipCum] + fabs(slip - slipPrev);
+
+ if (slipCum < properties[p_d0]) {
// if/else linear slip-weakening form of mu_f
mu_f = properties[p_coefS] -
(properties[p_coefS] - properties[p_coefD]) *
- stateVars[s_slipCum] / properties[p_d0];
+ slipCum / properties[p_d0];
} else {
mu_f = properties[p_coefD];
} // if/else
@@ -308,7 +320,7 @@
assert(_SlipWeakening::numStateVars == numStateVars);
const PylithScalar tolerance = 1.0e-12;
- if (slipRate > tolerance) {
+ if (slipRate > tolerance && !_forceHealing) {
const PylithScalar slipPrev = stateVars[s_slipPrev];
stateVars[s_slipPrev] = stateVars[s_slipCum];
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh 2012-01-27 01:09:27 UTC (rev 19485)
@@ -47,6 +47,12 @@
/// Destructor.
~SlipWeakening(void);
+ /** Compute properties from values in spatial database.
+ *
+ * @param flag True if forcing healing, false otherwise.
+ */
+ void forceHealing(const bool flag);
+
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -166,6 +172,8 @@
static const int db_slipCum;
static const int db_slipPrev;
+ bool _forceHealing; ///< Force healing.
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -146,7 +146,9 @@
CHECK_PETSC_ERROR(err);
// Account for censored cells
- const int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
+ int cellDepth = sieveMesh->depth();
+ err = MPI_Allreduce(&cellDepth, &cellDepth, 1, MPI_INT, MPI_MAX,
+ sieveMesh->comm());CHECK_PETSC_ERROR(err);
const int depth = (0 == label) ? cellDepth : labelId;
const std::string labelName = (0 == label) ?
((sieveMesh->hasLabel("censored depth")) ?
@@ -221,12 +223,12 @@
} catch (const std::exception& err) {
std::ostringstream msg;
msg << "Error while opening HDF5 file "
- << _filename << ".\n" << err.what();
+ << _hdf5Filename() << ".\n" << err.what();
throw std::runtime_error(msg.str());
} catch (...) {
std::ostringstream msg;
msg << "Unknown error while opening HDF5 file "
- << _filename << ".\n";
+ << _hdf5Filename() << ".\n";
throw std::runtime_error(msg.str());
} // try/catch
} // open
@@ -346,12 +348,12 @@
} catch (const std::exception& err) {
std::ostringstream msg;
msg << "Error while writing field '" << field.label() << "' at time "
- << t << " to HDF5 file '" << _filename << "'.\n" << err.what();
+ << t << " to HDF5 file '" << _hdf5Filename() << "'.\n" << err.what();
throw std::runtime_error(msg.str());
} catch (...) {
std::ostringstream msg;
msg << "Error while writing field '" << field.label() << "' at time "
- << t << " to HDF5 file '" << _filename << "'.\n";
+ << t << " to HDF5 file '" << _hdf5Filename() << "'.\n";
throw std::runtime_error(msg.str());
} // try/catch
} // writeVertexField
@@ -372,11 +374,14 @@
try {
const char* context = DataWriter<mesh_type, field_type>::_context.c_str();
+ PetscErrorCode err = 0;
const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
field.mesh().sieveMesh();
assert(!sieveMesh.isNull());
- const int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
+ int cellDepth = sieveMesh->depth();
+ err = MPI_Allreduce(&cellDepth, &cellDepth, 1, MPI_INT, MPI_MAX,
+ sieveMesh->comm());CHECK_PETSC_ERROR(err);
const int depth = (0 == label) ? cellDepth : labelId;
const std::string labelName = (0 == label) ?
((sieveMesh->hasLabel("censored depth")) ?
@@ -400,7 +405,6 @@
field.mesh().comm());
assert(fiberDim > 0);
- PetscErrorCode err = 0;
if (_timesteps.find(field.label()) == _timesteps.end())
_timesteps[field.label()] = 0;
@@ -434,12 +438,12 @@
} catch (const std::exception& err) {
std::ostringstream msg;
msg << "Error while writing field '" << field.label() << "' at time "
- << t << " to HDF5 file '" << _filename << "'.\n" << err.what();
+ << t << " to HDF5 file '" << _hdf5Filename() << "'.\n" << err.what();
throw std::runtime_error(msg.str());
} catch (...) {
std::ostringstream msg;
msg << "Error while writing field '" << field.label() << "' at time "
- << t << " to HDF5 file '" << _filename << "'.\n";
+ << t << " to HDF5 file '" << _hdf5Filename() << "'.\n";
throw std::runtime_error(msg.str());
} // try/catch
} // writeCellField
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -162,7 +162,9 @@
// Write cells
// Account for censored cells
- const int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
+ int cellDepth = sieveMesh->depth();
+ err = MPI_Allreduce(&cellDepth, &cellDepth, 1, MPI_INT, MPI_MAX,
+ sieveMesh->comm());CHECK_PETSC_ERROR(err);
const int depth = (0 == label) ? cellDepth : labelId;
const std::string labelName = (0 == label) ?
((sieveMesh->hasLabel("censored depth")) ?
@@ -436,14 +438,17 @@
try {
const char* context = DataWriter<mesh_type, field_type>::_context.c_str();
+ PetscErrorCode err = 0;
int rank = 0;
- MPI_Comm_rank(field.mesh().comm(), &rank);
+ err = MPI_Comm_rank(field.mesh().comm(), &rank);CHECK_PETSC_ERROR(err);
const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
field.mesh().sieveMesh();
assert(!sieveMesh.isNull());
- const int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
+ int cellDepth = sieveMesh->depth();
+ err = MPI_Allreduce(&cellDepth, &cellDepth, 1, MPI_INT, MPI_MAX,
+ sieveMesh->comm());CHECK_PETSC_ERROR(err);
const int depth = (0 == label) ? cellDepth : labelId;
const std::string labelName = (0 == label) ?
((sieveMesh->hasLabel("censored depth")) ?
@@ -458,7 +463,6 @@
assert(vector);
PetscViewer binaryViewer;
- PetscErrorCode err = 0;
const hid_t scalartype = (sizeof(double) == sizeof(PylithScalar)) ?
H5T_IEEE_F64BE : H5T_IEEE_F32BE;
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -133,15 +133,16 @@
err = VTKViewer::writeVertices(sieveMesh, _viewer);
CHECK_PETSC_ERROR(err);
//std::cout << "Wrote vertices for " << filename << std::endl;
- if (0 == label) {
- err = VTKViewer::writeElements(sieveMesh, _viewer);
- CHECK_PETSC_ERROR(err);
- } else {
- const std::string labelName =
- (sieveMesh->hasLabel("censored depth")) ? "censored depth" : "depth";
- err = VTKViewer::writeElements(sieveMesh, label, labelId, labelName, 0, _viewer);
- CHECK_PETSC_ERROR(err);
- } // if
+
+ const int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
+ const int depth = (!label) ? cellDepth : labelId;
+ const std::string cLabelName = (!label) ?
+ ((sieveMesh->hasLabel("censored depth")) ?
+ "censored depth" : "depth") : label;
+ const std::string vLabelName =
+ (sieveMesh->hasLabel("censored depth")) ? "censored depth" : "depth";
+ err = VTKViewer::writeElements(sieveMesh, cLabelName, depth, vLabelName, 0,
+ _viewer); CHECK_PETSC_ERROR(err);
//std::cout << "Wrote elements for " << filename << std::endl;
_wroteVertexHeader = false;
@@ -149,12 +150,12 @@
} catch (const std::exception& err) {
std::ostringstream msg;
msg << "Error while preparing for writing data to VTK file "
- << _filename << " at time " << t << ".\n" << err.what();
+ << _vtkFilename(t) << " at time " << t << ".\n" << err.what();
throw std::runtime_error(msg.str());
} catch (...) {
std::ostringstream msg;
msg << "Unknown error while preparing for writing data to VTK file "
- << _filename << " at time " << t << ".\n";
+ << _vtkFilename(t) << " at time " << t << ".\n";
throw std::runtime_error(msg.str());
} // try/catch
} // openTimeStep
@@ -222,12 +223,12 @@
} catch (const std::exception& err) {
std::ostringstream msg;
msg << "Error while writing field '" << field.label() << "' at time "
- << t << " to VTK file '" << _filename << "'.\n" << err.what();
+ << t << " to VTK file '" << _vtkFilename(t) << "'.\n" << err.what();
throw std::runtime_error(msg.str());
} catch (...) {
std::ostringstream msg;
msg << "Error while writing field '" << field.label() << "' at time "
- << t << " to VTK file '" << _filename << "'.\n";
+ << t << " to VTK file '" << _vtkFilename(t) << "'.\n";
throw std::runtime_error(msg.str());
} // try/catch
} // writeVertexField
@@ -246,26 +247,29 @@
typedef typename field_type::Mesh::RealSection RealSection;
try {
+ PetscErrorCode err = 0;
int rank = 0;
- MPI_Comm_rank(field.mesh().comm(), &rank);
+ err = MPI_Comm_rank(field.mesh().comm(), &rank);CHECK_PETSC_ERROR(err);
// Correctly handle boundary and fault meshes
// Cannot just use mesh->depth() because boundaries report the wrong thing
const ALE::Obj<SieveMesh>& sieveMesh = field.mesh().sieveMesh();
assert(!sieveMesh.isNull());
- const int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
- const int depth = (0 == label) ? cellDepth : labelId;
- const std::string labelName = (0 == label) ?
+ int cellDepth = sieveMesh->depth();
+ err = MPI_Allreduce(&cellDepth, &cellDepth, 1, MPI_INT, MPI_MAX,
+ sieveMesh->comm());CHECK_PETSC_ERROR(err);
+ const int depth = (!label) ? cellDepth : labelId;
+ const std::string labelName = (!label) ?
((sieveMesh->hasLabel("censored depth")) ?
"censored depth" : "depth") : label;
assert(!sieveMesh->getFactory().isNull());
const ALE::Obj<typename SieveMesh::numbering_type>& numbering =
sieveMesh->getFactory()->getNumbering(sieveMesh, labelName, depth);
assert(!numbering.isNull());
- assert(!sieveMesh->getLabelStratum(labelName, depth).isNull());
const ALE::Obj<RealSection>& section = field.section();
assert(!section.isNull());
+ assert(!sieveMesh->getLabelStratum(labelName, depth).isNull());
int fiberDimLocal =
(sieveMesh->getLabelStratum(labelName, depth)->size() > 0) ?
section->getFiberDimension(*sieveMesh->getLabelStratum(labelName, depth)->begin()) : 0;
@@ -273,10 +277,10 @@
MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX,
field.mesh().comm());
assert(fiberDim > 0);
+
const int enforceDim =
(field.vectorFieldType() != topology::FieldBase::VECTOR) ? fiberDim : 3;
- PetscErrorCode err = 0;
if (!_wroteCellHeader) {
err = PetscViewerASCIIPrintf(_viewer, "CELL_DATA %d\n",
numbering->getGlobalSize());
@@ -289,12 +293,12 @@
} catch (const std::exception& err) {
std::ostringstream msg;
msg << "Error while writing field '" << field.label() << "' at time "
- << t << " to VTK file '" << _filename << "'.\n" << err.what();
+ << t << " to VTK file '" << _vtkFilename(t) << "'.\n" << err.what();
throw std::runtime_error(msg.str());
} catch (...) {
std::ostringstream msg;
msg << "Error while writing field '" << field.label() << "' at time "
- << t << " to VTK file '" << _filename << "'.\n";
+ << t << " to VTK file '" << _vtkFilename(t) << "'.\n";
throw std::runtime_error(msg.str());
} // try/catch
} // writeCellField
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Solver.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -49,6 +49,7 @@
PetscFunctionBegin;
ierr = MatShellGetContext(mat, (void **) &ctx);CHKERRQ(ierr);
ierr = PCFieldSplitGetIS(ctx->pc, ctx->faultFieldName, &faultIS);CHKERRQ(ierr);
+ assert(faultIS);
ierr = ISEqual(isrow, faultIS, &isFaultRow);CHKERRQ(ierr);
ierr = ISEqual(iscol, faultIS, &isFaultCol);CHKERRQ(ierr);
if (isFaultRow && isFaultCol) {
@@ -210,7 +211,21 @@
_ctx.pc = *pc;
_ctx.A = jacobian.matrix();
- _ctx.faultFieldName = "3";
+ switch (spaceDim) {
+ case 1 :
+ _ctx.faultFieldName = "1";
+ break;
+ case 2 :
+ _ctx.faultFieldName = "2";
+ break;
+ case 3 :
+ _ctx.faultFieldName = "3";
+ break;
+ default:
+ assert(0);
+ throw std::logic_error("Unknown space dimension in "
+ "Problems::_setupFieldSplit().");
+ } // switch
_ctx.faultA = _jacobianPCFault;
} // if
} // _setupFieldSplit
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -98,6 +98,7 @@
err = SNESSetFromOptions(_snes); CHECK_PETSC_ERROR(err);
err = SNESLineSearchSet(_snes, lineSearch, (void*) formulation); CHECK_PETSC_ERROR(err);
+ err = SNESSetComputeInitialGuess(_snes, initialGuess, (void*) formulation); CHECK_PETSC_ERROR(err);
if (formulation->splitFields()) {
PetscKSP ksp = 0;
@@ -118,43 +119,6 @@
{ // solve
assert(0 != solution);
-#if 0
- // Update KSP operators with custom preconditioner if necessary.
- const ALE::Obj<RealSection>& solutionSection = solution->section();
- assert(!solutionSection.isNull());
- const ALE::Obj<SieveMesh>& sieveMesh = solution->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- if (solutionSection->getNumSpaces() > sieveMesh->getDimension() &&
- 0 != _jacobianPCFault) {
- PetscKSP ksp = 0;
- PetscPC pc = 0;
- PetscKSP *ksps = 0;
- PetscMat A = 0;
- PylithInt num = 0;
-
- PetscErrorCode err = 0;
- err = SNESGetKSP(_snes, &ksp); CHECK_PETSC_ERROR(err);
- err = KSPSetUp(ksp); CHECK_PETSC_ERROR(err);
- err = KSPGetPC(ksp, &pc); CHECK_PETSC_ERROR(err);
- err = PCFieldSplitGetSubKSP(pc, &num, &ksps); CHECK_PETSC_ERROR(err);
- assert(solutionSection->getNumSpaces() == num);
-
-#if 0 // debugging
- std::cout << "Preconditioner Matrix" << std::endl;
- MatView(_jacobianPCFault, PETSC_VIEWER_STDOUT_WORLD);
-#endif
-
-
- MatStructure flag;
- err = KSPGetOperators(ksps[num-1], &A,
- PETSC_NULL, &flag); CHECK_PETSC_ERROR(err);
- err = PetscObjectReference((PetscObject) A); CHECK_PETSC_ERROR(err);
- err = KSPSetOperators(ksps[num-1], A, _jacobianPCFault,
- flag); CHECK_PETSC_ERROR(err);
- err = PetscFree(ksps); CHECK_PETSC_ERROR(err);
- } // if
-#endif
-
const int solveEvent = _logger->eventId("SoNl solve");
const int scatterEvent = _logger->eventId("SoNl scatter");
_logger->eventBegin(solveEvent);
@@ -361,11 +325,16 @@
ierr = PetscViewerASCIIPrintf(snes->ls_monitor," Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",(double)fnorm,(double)*gnorm,(double)*ynorm,(double)minlambda,(double)lambda,(double)initslope);CHKERRQ(ierr);
ierr = PetscViewerASCIISubtractTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
}
- *flag = PETSC_FALSE; // DIVERGED_LINE_SEARCH
assert(lsctx);
Formulation* formulation = (Formulation*) lsctx;
assert(formulation);
formulation->printState(&w, &g, &x, &y);
+#if 0 // Original PETSc code
+ *flag = PETSC_FALSE; // DIVERGED_LINE_SEARCH
+#else
+ std::cerr << "WARNING: Line search diverged ... continuing nonlinear iterations anyway in hopes that solution will converge anyway."
+ << std::endl;
+#endif
break;
}
t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope;
@@ -377,24 +346,7 @@
if (a == 0.0) {
lambdatemp = -initslope/(2.0*b);
} else {
-
- // MATT: Check this
- //
- // Temporary fix by Brad to keep lambda in proper
- // range. Necessary due to underflow and overflow of a, b, c,
- // and d.
-#if 0
lambdatemp = (-b + PetscSqrtScalar(d))/(3.0*a);
-#else
- if ((-b + PetscSqrtScalar(d) > 0.0 && a > 0.0) ||
- (-b + PetscSqrtScalar(d) < 0.0 && a < 0.0)) {
- lambdatemp = (-b + PetscSqrtScalar(d))/(3.0*a);
- } else {
- lambdatemp = 0.05*lambda;
- } // else
-#endif
-
-
} // if/else
lambdaprev = lambda;
gnormprev = *gnorm;
@@ -405,8 +357,6 @@
else
lambda = lambdatemp;
-
-
ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr);
if (snes->nfuncs >= snes->max_funcs) {
ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr);
@@ -457,12 +407,6 @@
ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
// ======================================================================
- // Code to constrain solution space.
- assert(lsctx);
- Formulation* formulation = (Formulation*) lsctx;
- assert(formulation);
- formulation->constrainSolnSpace(&w);
-
ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
if (snes->domainerror) {
ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr);
@@ -482,6 +426,18 @@
} // lineSearch
// ----------------------------------------------------------------------
+// Generic C interface for customized PETSc initial guess.
+#undef __FUNCT__
+#define __FUNCT__ "initialGuess"
+PetscErrorCode
+pylith::problems::SolverNonlinear::initialGuess(PetscSNES snes,
+ PetscVec initialGuessVec,
+ void *lsctx)
+{ // initialGuess
+ VecSet(initialGuessVec, 0.0);
+} // initialGuess
+
+// ----------------------------------------------------------------------
// Initialize logger.
void
pylith::problems::SolverNonlinear::_initializeLogger(void)
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.hh 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.hh 2012-01-27 01:09:27 UTC (rev 19485)
@@ -141,6 +141,18 @@
PetscReal *gnorm,
PetscBool *flag);
+ /** Generic C interface for customized PETSc initial guess.
+ *
+ * @param snes PETSc SNES solver.
+ * @param initialGuessVec PETSc vector for initial guess.
+ * @param lsctx Optional context for line search (not used here).
+ * @returns PETSc error code.
+ */
+ static
+ PetscErrorCode initialGuess(PetscSNES snes,
+ PetscVec initialGuessVec,
+ void *lsctx);
+
// PRIVATE METHODS //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/MeshRefiner.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/MeshRefiner.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/MeshRefiner.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -22,6 +22,8 @@
#include "MeshOrder.hh" // USES MeshOrder
+#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
+
#include <cassert> // USES assert()
// ----------------------------------------------------------------------
@@ -51,11 +53,34 @@
cellrefiner_type& refiner)
{ // refine
assert(!mesh.isNull());
+
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ //logger.setDebug(1);
+ logger.stagePush("RefinedMesh");
+
if (mesh->hasLabel("censored depth")) {
_refineCensored(newMesh, mesh, refiner);
} else {
_refine(newMesh, mesh, refiner);
} // if/else
+
+ logger.stagePush("RefinedMeshStratification");
+ _stratify(newMesh);
+ logger.stagePop();
+
+ logger.stagePush("RefinedMeshOverlap");
+ _calcNewOverlap(newMesh, mesh, refiner);
+ logger.stagePop();
+
+ logger.stagePush("RefinedMeshIntSections");
+ _createIntSections(newMesh, mesh, refiner);
+ logger.stagePop();
+
+ logger.stagePush("RefinedMeshLabels");
+ _createLabels(newMesh, mesh, refiner);
+ logger.stagePop();
+
+ logger.stagePop(); // RefinedMesh
} // refine
// ----------------------------------------------------------------------
@@ -73,7 +98,6 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
//logger.setDebug(1);
- logger.stagePush("RefinedMesh");
logger.stagePush("RefinedMeshCreation");
// Calculate order in old mesh.
@@ -195,27 +219,6 @@
refiner.setCoordsNewVertices(newCoordinates, coordinates);
logger.stagePop();
- logger.stagePush("RefinedMeshStratification");
-
- _stratify(newMesh);
-
- logger.stagePop();
- logger.stagePush("RefinedMeshOverlap");
-
- _calcNewOverlap(newMesh, mesh, refiner);
-
- logger.stagePop();
- logger.stagePush("RefinedMeshIntSections");
-
- _createIntSections(newMesh, mesh, refiner);
-
- logger.stagePop();
- logger.stagePush("RefinedMeshLabels");
-
- _createLabels(newMesh, mesh, refiner);
-
- logger.stagePop();
- logger.stagePop(); // Mesh
} // _refine
// ----------------------------------------------------------------------
@@ -432,11 +435,6 @@
refiner.setCoordsNewVertices(newCoordinates, coordinates);
- _stratify(newMesh);
- _calcNewOverlap(newMesh, mesh, refiner);
- _createIntSections(newMesh, mesh, refiner);
- _createLabels(newMesh, mesh, refiner);
-
// Create sensored depth
const ALE::Obj<SieveFlexMesh::label_type>& censoredLabel = newMesh->createLabel("censored depth");
assert(!censoredLabel.isNull());
@@ -709,10 +707,12 @@
if (localPoint >= oldVerticesStartNormal && localPoint < oldVerticesStartFault) {
assert(remotePoint >= oldVerticesStartNormalP[rank] && remotePoint < oldVerticesStartFaultP[rank]);
+ assert(remotePoint+remoteOffsetNormal >= 0);
newSendOverlap->addArrow(localPoint+localOffsetNormal, rank, remotePoint+remoteOffsetNormal);
} else {
assert(localPoint >= oldVerticesStartFault);
assert(remotePoint >= oldVerticesStartFaultP[rank]);
+ assert(remotePoint+remoteOffsetFault >= 0);
newSendOverlap->addArrow(localPoint+localOffsetFault, rank, remotePoint+remoteOffsetFault);
}
} // for
@@ -756,10 +756,12 @@
if (localPoint >= oldVerticesStartNormal && localPoint < oldVerticesStartFault) {
assert(remotePoint >= oldVerticesStartNormalP[rank] && remotePoint < oldVerticesStartFaultP[rank]);
+ assert(remotePoint+remoteOffsetNormal >= 0);
newSendOverlap->addArrow(localPoint+localOffsetNormal, rank, remotePoint+remoteOffsetNormal);
} else {
assert(localPoint >= oldVerticesStartFault);
assert(remotePoint >= oldVerticesStartFaultP[rank]);
+ assert(remotePoint+remoteOffsetFault >= 0);
newSendOverlap->addArrow(localPoint+localOffsetFault, rank, remotePoint+remoteOffsetFault);
}
}
@@ -799,12 +801,79 @@
// We have to do flexible assembly since we add the new vertices separately
newSendOverlap->assemble();
newRecvOverlap->assemble();
+
+ // Verify size of new send/recv overlaps are at least as big as the
+ // original ones.
+ PetscErrorCode err = 0;
+ if (newSendOverlap->getNumRanks() != sendOverlap->getNumRanks() ||
+ newRecvOverlap->getNumRanks() != recvOverlap->getNumRanks() ||
+ newSendOverlap->getNumRanks() != newRecvOverlap->getNumRanks()) {
+
+ throw std::logic_error("Error in constructing new overlaps during mesh "
+ "refinement.\nMismatch in number of ranks.");
+ } // if
+
+ const int numRanks = newSendOverlap->getNumRanks();
+ for (int isend=0; isend < numRanks; ++isend) {
+ int rank = 0;
+ err = newSendOverlap->getRank(isend, &rank); CHECK_PETSC_ERROR(err);
+ int irecv = 0;
+ err = newRecvOverlap->getRankIndex(rank, &irecv);
+
+ int sendRank = 0;
+ int recvRank = 0;
+ err = sendOverlap->getRank(isend, &sendRank); CHECK_PETSC_ERROR(err);
+ err = recvOverlap->getRank(irecv, &recvRank); CHECK_PETSC_ERROR(err);
+ if (rank != sendRank || rank != recvRank) {
+ throw std::logic_error("Error in constructing new overlaps during mesh "
+ "refinement.\nMismatch in ranks.");
+ } // if
+
+ int newSendNumPoints = 0;
+ int oldSendNumPoints = 0;
+ int newRecvNumPoints = 0;
+ int oldRecvNumPoints = 0;
+ err = sendOverlap->getNumPointsByRank(rank, &oldSendNumPoints); CHECK_PETSC_ERROR(err);
+ err = newSendOverlap->getNumPointsByRank(rank, &newSendNumPoints); CHECK_PETSC_ERROR(err);
+ err = recvOverlap->getNumPointsByRank(rank, &oldRecvNumPoints); CHECK_PETSC_ERROR(err);
+ err = newRecvOverlap->getNumPointsByRank(rank, &newRecvNumPoints); CHECK_PETSC_ERROR(err);
+ if (newSendNumPoints < oldSendNumPoints || newRecvNumPoints < oldRecvNumPoints) {
+ throw std::logic_error("Error in constructing new overlaps during mesh "
+ "refinement.\nInvalid size for new overlaps.");
+ } // if
+ } // for
+
if (mesh->debug()) {
sendOverlap->view("OLD SEND OVERLAP");
recvOverlap->view("OLD RECV OVERLAP");
newSendOverlap->view("NEW SEND OVERLAP");
newRecvOverlap->view("NEW RECV OVERLAP");
- }
+
+ int rank = 0;
+ int nprocs = 0;
+ MPI_Comm_rank(mesh->comm(), &rank);
+ MPI_Comm_size(mesh->comm(), &nprocs);
+ MPI_Barrier(mesh->comm());
+ for (int i=0; i < nprocs; ++i) {
+ if (rank == i) {
+ int numNeighbors = newSendOverlap->getNumRanks();
+ assert(numNeighbors == newRecvOverlap->getNumRanks());
+ for (int j=0; j < numNeighbors; ++j) {
+ int sendRank = 0;
+ int recvRank = 0;
+ err = newSendOverlap->getRank(j, &sendRank);CHECK_PETSC_ERROR(err);
+ err = newRecvOverlap->getRank(j, &recvRank);CHECK_PETSC_ERROR(err);
+ std::cout << "["<<rank<<"]: "
+ << "send: " << sendRank
+ << ", npts: " << newSendOverlap->getNumPoints(j)
+ << "; recv: " << recvRank
+ << ", npts: " << newRecvOverlap->getNumPoints(j)
+ << std::endl;
+ } // for
+ } // if
+ MPI_Barrier(mesh->comm());
+ } // for
+ } // if
} // _calcNewOverlap
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/RefineEdges2.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/RefineEdges2.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/RefineEdges2.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -197,7 +197,7 @@
std::set_intersection(leftRanks.begin(), leftRanks.end(), rightRanks.begin(), rightRanks.end(),
std::insert_iterator<std::set<int> >(ranks, ranks.begin()));
-#if 0
+#if 0 // DEBUGGING
std::cout << "[" << myrank << "] Checking edge " << e_iter->first << std::endl;
for(std::set<int>::const_iterator r_iter = leftRanks.begin(); r_iter != leftRanks.end(); ++r_iter) {
std::cout << "[" << myrank << "] left rank " << *r_iter << std::endl;
@@ -212,7 +212,11 @@
newVerticesSection->addFiberDimension(edgeMin+localMinOffset, 1);
for(std::set<int>::const_iterator r_iter = ranks.begin(); r_iter != ranks.end(); ++r_iter) {
bndryEdgeToRank[e_iter->first].push_back(*r_iter);
- //std::cout << "[" << myrank << "] Added edge " << e_iter->first << " with rank " << *r_iter << std::endl;
+#if 0 // DEBUGGING
+ const point_type edgeMax = std::max(e_iter->first.first, e_iter->first.second);
+ const int localMaxOffset = (orderOldMesh.verticesNormal().hasPoint(edgeMax)) ? localNormalOffset : localCensoredOffset;
+ std::cout << "[" << myrank << "] Added edge " << e_iter->first << " now (" << edgeMin+localMinOffset << ", " << edgeMax+localMaxOffset << ") with rank " << *r_iter << std::endl;
+#endif
} // for
} // if
} // if
@@ -245,8 +249,10 @@
Obj<ALE::Section<overlap_point_type, EdgeType> > overlapVertices = new ALE::Section<overlap_point_type, EdgeType>(oldMesh->comm());
ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newVerticesSection, overlapVertices);
- //newVerticesSection->view("NEW VERTICES");
- //overlapVertices->view("OVERLAP VERTICES");
+#if 0 // DEBUGGING
+ newVerticesSection->view("NEW VERTICES");
+ overlapVertices->view("OVERLAP VERTICES");
+#endif
// Merge by translating edge to local points, finding edge in _edgeToVertex, and adding (local new vetex, remote new vertex) to overlap
for(std::map<EdgeType, std::vector<int> >::const_iterator e_iter = bndryEdgeToRank.begin(); e_iter != bndryEdgeToRank.end(); ++e_iter) {
@@ -284,9 +290,10 @@
break;
} // if
} // for
-#if 0
+#if 0 // DEBUGGING
if (-1 == newRemotePoint) {
- std::cout << "remoteLeft: " << remoteLeft
+ std::cout << "["<< myrank << "] DISMISSING newLocalPoint: " << newLocalPoint
+ << ", remoteLeft: " << remoteLeft
<< ", remoteRight: " << remoteRight
<< ", rank: " << rank
<< ", remoteSize: " << remoteSize
@@ -300,13 +307,6 @@
} // if
} // for
} // for
-
-#if 0
- oldSendOverlap->view("OLD SEND OVERLAP");
- oldRecvOverlap->view("OLD RECV OVERLAP");
- newSendOverlap->view("NEW SEND OVERLAP");
- newRecvOverlap->view("NEW RECV OVERLAP");
-#endif
} // overlapAddNewVertces
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/RefineFace4Edges2.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/RefineFace4Edges2.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/RefineFace4Edges2.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -345,13 +345,6 @@
}
} // for
} // for
-
-#if 0
- oldSendOverlap->view("OLD SEND OVERLAP");
- oldRecvOverlap->view("OLD RECV OVERLAP");
- newSendOverlap->view("NEW SEND OVERLAP");
- newRecvOverlap->view("NEW RECV OVERLAP");
-#endif
} // overlapAddNewVertces
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/RefineUniform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/RefineUniform.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/RefineUniform.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -119,6 +119,7 @@
break;
default :
+ assert(0);
throw std::logic_error("Unknown dimension.");
} // switch
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/RefineVol8Face4Edges2.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/RefineVol8Face4Edges2.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/RefineVol8Face4Edges2.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -588,8 +588,7 @@
} // if
} // for
assert(localPoint >= orderNewMesh.verticesNormal().min() && localPoint < orderNewMesh.verticesNormal().max());
-#if 0
- // Debugging for fault edge problem
+#if 0 // Debugging for fault edge problem
assert(remotePoint >= 0);
#endif
if (remotePoint >= 0) {
@@ -598,13 +597,6 @@
}
} // for
} // for
-
-#if 0 // debuggin
- oldSendOverlap->view("OLD SEND OVERLAP");
- oldRecvOverlap->view("OLD RECV OVERLAP");
- newSendOverlap->view("NEW SEND OVERLAP");
- newRecvOverlap->view("NEW RECV OVERLAP");
-#endif
} // overlapAddNewVertces
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/ReverseCuthillMcKee.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/ReverseCuthillMcKee.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/ReverseCuthillMcKee.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -33,7 +33,7 @@
void
pylith::topology::ReverseCuthillMcKee::reorder(topology::Mesh* mesh)
{ // reorder
- assert(0 != mesh);
+ assert(mesh);
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
//logger.setDebug(1);
Modified: short/3D/PyLith/trunk/modulesrc/friction/SlipWeakening.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/friction/SlipWeakening.i 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/modulesrc/friction/SlipWeakening.i 2012-01-27 01:09:27 UTC (rev 19485)
@@ -36,6 +36,12 @@
/// Destructor.
~SlipWeakening(void);
+ /** Compute properties from values in spatial database.
+ *
+ * @param flag True if forcing healing, false otherwise.
+ */
+ void forceHealing(const bool flag);
+
// PROTECTED METHODS //////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/pylith/apps/PyLithApp.py
===================================================================
--- short/3D/PyLith/trunk/pylith/apps/PyLithApp.py 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/pylith/apps/PyLithApp.py 2012-01-27 01:09:27 UTC (rev 19485)
@@ -63,7 +63,14 @@
factory=MemoryLogger)
perfLogger.meta['tip'] = "Performance and memory logging."
+
+ typos = pyre.inventory.str("typos", default="pedantic",
+ validator=pyre.inventory.choice(['relaxed', 'strict', 'pedantic']))
+ typos.meta['tip'] = "Specifies the handling of unknown properties and " \
+ "facilities"
+
+
# PUBLIC METHODS /////////////////////////////////////////////////////
def __init__(self, name="pylithapp"):
Modified: short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/Integrator.py 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/pylith/feassemble/Integrator.py 2012-01-27 01:09:27 UTC (rev 19485)
@@ -123,6 +123,7 @@
"""
return
+
# PRIVATE METHODS ////////////////////////////////////////////////////
def _setupLogging(self):
Modified: short/3D/PyLith/trunk/pylith/friction/RateStateAgeing.py
===================================================================
--- short/3D/PyLith/trunk/pylith/friction/RateStateAgeing.py 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/pylith/friction/RateStateAgeing.py 2012-01-27 01:09:27 UTC (rev 19485)
@@ -76,7 +76,7 @@
'cell': \
{'info': [],
'data': []}}
- self._loggingPrefix = "FrRSA "
+ self._loggingPrefix = "FrRSAg "
return
Modified: short/3D/PyLith/trunk/pylith/friction/SlipWeakening.py
===================================================================
--- short/3D/PyLith/trunk/pylith/friction/SlipWeakening.py 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/pylith/friction/SlipWeakening.py 2012-01-27 01:09:27 UTC (rev 19485)
@@ -33,6 +33,28 @@
Factory: friction_model.
"""
+ # INVENTORY //////////////////////////////////////////////////////////
+
+ class Inventory(FrictionModel.Inventory):
+ """
+ Python object for managing SlipWeakening facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing FrictionModel facilities and properties.
+ ##
+ ## \b Properties
+ ## @li \b force_healing Force healing after every time step.
+ ##
+ ## \b Facilities
+ ## @li None
+
+ import pyre.inventory
+
+ forceHealing = pyre.inventory.bool("force_healing", default=False)
+ forceHealing.meta['tip'] = "Force healing after every time step."
+
+
# PUBLIC METHODS /////////////////////////////////////////////////////
def __init__(self, name="slipweakening"):
@@ -51,12 +73,25 @@
'cell': \
{'info': [],
'data': []}}
- self._loggingPrefix = "FrSW "
+ self._loggingPrefix = "FrSlWk "
return
# PRIVATE METHODS ////////////////////////////////////////////////////
+ def _configure(self):
+ """
+ Setup members using inventory.
+ """
+ try:
+ FrictionModel._configure(self)
+ ModuleSlipWeakening.forceHealing(self, self.inventory.forceHealing)
+ except ValueError, err:
+ aliases = ", ".join(self.aliases)
+ raise ValueError("Error while configuring friction model "
+ "(%s):\n%s" % (aliases, err.message))
+ return
+
def _createModuleObj(self):
"""
Call constructor for module object for access to C++ object.
Modified: short/3D/PyLith/trunk/pylith/friction/TimeWeakening.py
===================================================================
--- short/3D/PyLith/trunk/pylith/friction/TimeWeakening.py 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/pylith/friction/TimeWeakening.py 2012-01-27 01:09:27 UTC (rev 19485)
@@ -50,7 +50,7 @@
'cell': \
{'info': [],
'data': []}}
- self._loggingPrefix = "FrStat "
+ self._loggingPrefix = "FrTmWk "
return
Modified: short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg 2012-01-27 01:09:27 UTC (rev 19485)
@@ -47,7 +47,7 @@
[pylithapp.timedependent.implicit.time_step]
total_time = 14.0*s
-dt = 0.05*s
+dt = 0.1*s
# ----------------------------------------------------------------------
# materials
@@ -98,8 +98,8 @@
# KSP
ksp_rtol = 1.0e-14
-ksp_atol = 1.0e-18
-ksp_max_it = 500
+ksp_atol = 1.0e-15
+ksp_max_it = 50
ksp_gmres_restart = 20
ksp_monitor = true
@@ -107,9 +107,9 @@
ksp_converged_reason = true
# SNES
-snes_rtol = 1.0e-12
-snes_atol = 1.0e-16
-snes_max_it = 500
+snes_rtol = 1.0e-14
+snes_atol = 1.0e-12
+snes_max_it = 100
snes_monitor = true
snes_ls_monitor = true
Modified: short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate.cfg 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate.cfg 2012-01-27 01:09:27 UTC (rev 19485)
@@ -46,8 +46,6 @@
db_initial_tractions.values = [traction-shear,traction-normal]
db_initial_tractions.data = [0.0*MPa, -10.0*MPa]
-zero_tolerance = 1.0e-14
-
# ----------------------------------------------------------------------
# PETSc
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate_stable.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate_stable.cfg 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate_stable.cfg 2012-01-27 01:09:27 UTC (rev 19485)
@@ -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/trunk/tests/2d/frictionslide/ratestate_weak.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate_weak.cfg 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate_weak.cfg 2012-01-27 01:09:27 UTC (rev 19485)
@@ -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-9
# Set the friction model parameters.
# reference coefficient of friction: 0.6
@@ -31,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/trunk/tests/2d/frictionslide/tension.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/frictionslide/tension.cfg 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/tests/2d/frictionslide/tension.cfg 2012-01-27 01:09:27 UTC (rev 19485)
@@ -88,13 +88,14 @@
[pylithapp.timedependent.interfaces.fault]
friction = pylith.friction.SlipWeakening
friction.label = Slip weakening
+friction.force_healing = True
# Set slip-weakening friction model parameters using a uniform DB. Set the
# parameters as follows:
friction.db_properties = spatialdata.spatialdb.UniformDB
friction.db_properties.label = Slip weakening
friction.db_properties.values = [static-coefficient,dynamic-coefficient,slip-weakening-parameter,cohesion]
-friction.db_properties.data = [0.6,0.0,1.0*mm,0.0*Pa]
+friction.db_properties.data = [0.6,0.2,1.0*mm,0.0*Pa]
db_initial_tractions = spatialdata.spatialdb.UniformDB
db_initial_tractions.label = Initial fault tractions
Modified: short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg 2012-01-27 01:09:27 UTC (rev 19485)
@@ -74,7 +74,7 @@
db_change = spatialdata.spatialdb.UniformDB
db_change.label = Amplitude of Dirichlet BC on +x
db_change.values = [displacement-x,displacement-y,displacement-z,change-start-time]
-db_change.data = [-0.2*m,-0.4*m,0.0*m,0.0*hour]
+db_change.data = [-0.5*m,-0.2*m,0.0*m,0.0*hour]
th_change = spatialdata.spatialdb.TimeHistory
th_change.label = Time history of tidal load
@@ -89,7 +89,7 @@
db_change = spatialdata.spatialdb.UniformDB
db_change.label = Amplitude of Dirichlet BC on -x
db_change.values = [displacement-x,displacement-y,displacement-z,change-start-time]
-db_change.data = [+0.2*m,+0.4*m,0.0*m,0.0*hour]
+db_change.data = [+0.5*m,+0.2*m,0.0*m,0.0*hour]
th_change = spatialdata.spatialdb.TimeHistory
th_change.label = Time history of tidal load
@@ -114,10 +114,11 @@
[pylithapp.timedependent.interfaces.fault]
label = fault
-zero_tolerance = 1.0e-8
+zero_tolerance = 1.0e-10
friction = pylith.friction.SlipWeakening
friction.label = Slip weakening
+friction.force_healing = True
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 2
@@ -170,8 +171,8 @@
sub_pc_factor_shift_type = nonzero
# Convergence parameters.
-ksp_rtol = 1.0e-8
-ksp_atol = 1.0e-10
+ksp_rtol = 1.0e-10
+ksp_atol = 1.0e-12
ksp_max_it = 100
ksp_gmres_restart = 50
@@ -181,8 +182,8 @@
ksp_converged_reason = true
# Nonlinear solver monitoring options.
-snes_rtol = 1.0e-8
-snes_atol = 1.0e-7
+snes_rtol = 1.0e-10
+snes_atol = 1.0e-9
snes_max_it = 100
snes_monitor = true
snes_ls_monitor = true
Copied: short/3D/PyLith/trunk/tests/3d/refine (from rev 19484, short/3D/PyLith/branches/v1.6-stable/tests/3d/refine)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-01-27 00:23:13 UTC (rev 19484)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-01-27 01:09:27 UTC (rev 19485)
@@ -216,6 +216,7 @@
const PylithScalar dt = 0.01;
fault.timeStep(dt);
fault.constrainSolnSpace(&fields, t, jacobian);
+ fault.updateStateVars(t, &fields);
{ // Check solution values
// No change to Lagrange multipliers for stick case.
@@ -334,6 +335,7 @@
const PylithScalar dt = 0.01;
fault.timeStep(dt);
fault.constrainSolnSpace(&fields, t, jacobian);
+ fault.updateStateVars(t, &fields);
{ // Check solution values
// Lagrange multipliers should be adjusted according to friction
@@ -458,6 +460,7 @@
const PylithScalar dt = 0.01;
fault.timeStep(dt);
fault.constrainSolnSpace(&fields, t, jacobian);
+ fault.updateStateVars(t, &fields);
//residual.view("RESIDUAL"); // DEBUGGING
More information about the CIG-COMMITS
mailing list