[cig-commits] r15640 - short/3D/PyLith/branches/pylith-friction/playpen/friction

[brad at geodynamics.org]

Author: brad
Date: 2009-08-31 20:26:42 -0700 (Mon, 31 Aug 2009)
New Revision: 15640

Modified:
   short/3D/PyLith/branches/pylith-friction/playpen/friction/spring_example.py
   short/3D/PyLith/branches/pylith-friction/playpen/friction/twoquad4.cfg
Log:
Improved spring example. Added description.

Modified: short/3D/PyLith/branches/pylith-friction/playpen/friction/spring_example.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/playpen/friction/spring_example.py	2009-08-31 23:50:46 UTC (rev 15639)
+++ short/3D/PyLith/branches/pylith-friction/playpen/friction/spring_example.py	2009-09-01 03:26:42 UTC (rev 15640)
@@ -1,15 +1,38 @@
+# Simple spring system to test friction implementation. Sliding occurs
+# between DOF 2 and 3. We specify the displacement at DOF 5.
 #
+# |--\/\/\--*--\/\/\--*--\/\/\--| |--\/\/\--*--\/\/\--*
+#      k0  u0    k1  u1    k2  u2 u3   k3  u4    k4   u5
 #
+# Want to solve A u(t+dt) = b(t+dt).
 #
+# Use displacement increment formulation, u(t+dt) = u(t) + du(t).
 #
+# Solve for displacement increment using A du(t) = b(t+dt) - A u(t).
+#
+# Residual at time t+dt is
+# r(t+dt) = b(t+dt) - A (u(t) + du(t)).
+#
+# Constraint is the u3 - u2 = 0, unless force exceeds static friction
+# (fy). If force exceeds friction, the friction is equal fy.
+#
+# Use Lagrange multipliers to enforce constraint.
+#
+# A u = b, becomes
+#
+# [ A C^T ] [ u ] = [ b ]
+# [ C  0  ] [ l ]   [ d ]
+#
+# where C is the constraint matrix, l is the vector of Lagrange
+# multipliers, and d is the vector of slip.
 
 import numpy
 
-# Spring stiffness
-k = (1.0, 1.0, 1.0, 1.0, 1.0)
+# Spring stiffness [k0, k1, ..., kN]
+k = (1.0, 1.0, 2.0, 1.0, 1.0)
 
-# Prescribed displacement
-u0 = 2.5
+# Prescribed displacement, u5
+u5 = 2.5
 
 # Static friction force
 fy = 0.4
@@ -17,6 +40,8 @@
 # Jacobian matrix of system
 # [ K  C^T ]
 # [ C  0   ]
+# K is stiffness matrix
+# C is the constraint matrix 
 A = numpy.array([[k[0]+k[1], -k[1], 0, 0, 0, 0],
                  [-k[1], k[1]+k[2], -k[2], 0, 0, 0],
                  [0, -k[2], k[2], 0, 0, -1],
@@ -34,7 +59,7 @@
     Calculate residual
     """
     residual = b - numpy.dot(A, disp)
-    residual[4] += k[4]*u0
+    residual[4] += k[4]*u5
     return residual
 
 
@@ -43,17 +68,15 @@
     """
     Calculate friction.
     """
-    global slipPrev
-    global s
-    global sset
-    dispTpdt = disp+incr
-    lm = dispTpdt[5]
+    dispTpdt = disp + incr
+    lm = dispTpdt[5] # Lagrange multiplier
     if lm > fy or (lm < fy and b[5] > 0.0):
         incr[5] = fy - disp[5] # Adjust Lagrange multiplier in solution
-        slip = 2*(lm-fy)
-        b[5] += slip # Update slip estimate
-        incr[2] -= 0.5*slip
-        incr[3] += 0.5*slip
+        #slipIncr = 2*(lm-fy)
+        slipIncr = (k[2] + k[3])*(lm-fy) - (k[2]*b[2] - k[3]*b[3])
+        b[5] += slipIncr # Update slip estimate
+        incr[2] -= 0.5*slipIncr
+        incr[3] += 0.5*slipIncr
 
 
 # ----------------------------------------------------------------------
@@ -61,26 +84,25 @@
     """
     Solve for increment of displacment given jacobian and current disp.
     """
-    i = 0
-    while numpy.sum(numpy.abs(residual)) > 1.0e-6 and i < 40:
-        print "Interation: %d" % i
-        ii = numpy.dot(Ai, residual)
-        incr += ii
+    iter = 0
+    while numpy.sum(numpy.abs(residual)) > 1.0e-8 and iter < 40:
+        print "Interation: %d" % iter
+        dincr = numpy.dot(Ai, residual) # Increment to disp increment.
+        incr += dincr
         calcFriction(disp, incr)
         residual = reformResidual(disp+incr)
         print "Disp(t):",disp
         print "Incr(t):",incr
         print "b:",b
         print "Residual:",residual
-        print "ii:",ii
-        i += 1
+        print "dincr:",dincr
+        iter += 1
     return
 
 # ----------------------------------------------------------------------
 # main
 incr = 0*disp
-residual = reformResidual(disp+incr)
+residual = reformResidual(disp + incr)
 solve(incr, A, residual, disp)
 disp += incr
 print "Solution:",disp
-print A

Modified: short/3D/PyLith/branches/pylith-friction/playpen/friction/twoquad4.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-friction/playpen/friction/twoquad4.cfg	2009-08-31 23:50:46 UTC (rev 15639)
+++ short/3D/PyLith/branches/pylith-friction/playpen/friction/twoquad4.cfg	2009-09-01 03:26:42 UTC (rev 15640)
@@ -58,6 +58,7 @@
 dimension = 2
 normalizer.length_scale = 1.0*m
 formulation = pylith.problems.Implicit
+formulation.solver = pylith.problems.SolverNonlinear
 
 # Set bc to an array with 2 boundary conditions: 'x_neg' and 'x_pos'.
 bc = [x_neg,x_pos]
@@ -186,6 +187,11 @@
 ksp_gmres_restart = 50
 #start_in_debugger = true
 
+snes_monitor = true
+snes_view = true
+ksp_converged_reason = true
+snes_converged_reason = true
+
 # ----------------------------------------------------------------------
 # output
 # ----------------------------------------------------------------------