[cig-commits] r19977 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src

[surendra at geodynamics.org]

Author: surendra
Date: 2012-04-23 18:41:02 -0700 (Mon, 23 Apr 2012)
New Revision: 19977

Modified:
   seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Switched to formulation in terms of tractions for Newton-Raphson search.  Added comments that were removed by mistake in previous commit

Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90	2012-04-24 01:17:41 UTC (rev 19976)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90	2012-04-24 01:41:02 UTC (rev 19977)
@@ -325,6 +325,9 @@
   ! NOTE: same Bi on both sides, see note above
   allocate(bc%Z(bc%nglob))
   bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
+  ! WARNING: In non-split nodes at fault edges M is assembled across the fault.
+  ! hence invM1+invM2=2/(M1+M2) instead of 1/M1+1/M2
+  ! In a symmetric mesh (M1=M2) Z will be twice its intended value
 
   allocate(bc%T(3,bc%nglob))
   allocate(bc%D(3,bc%nglob))
@@ -614,6 +617,7 @@
 !---------------------------------------------------------------------
 ! adds a value to a fault parameter inside an area with prescribed shape
 subroutine init_2d_distribution(a,coord,iin,n)
+!JPA refactor: background value shuld be an argument
 
   real(kind=CUSTOM_REAL), intent(inout) :: a(:)
   real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
@@ -695,6 +699,9 @@
 !=====================================================================
 ! adds boundary term Bt into Force array for each fault.
 !
+! NOTE: On non-split nodes at fault edges, dD=dV=dA=0
+! and the net contribution of B*T is =0
+!
 subroutine bc_dynflt_set3d_all(F,Vel,Dis)
 
   real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: Vel,Dis
@@ -806,6 +813,7 @@
     TxExt = DTau0 * bc%asp%Fload * GLoad
     T(1,:) = T(1,:) + TxExt
 
+    !JPA the solver below can be refactored into a loop with two passes
     Vf = Vnorm_old
     theta_old = bc%rsf%theta
     call rsf_update_state(Vf,bc%dt,bc%rsf)
@@ -910,6 +918,9 @@
   da(2,:) = bc%invM2*f(2,bc%ibulk2)-bc%invM1*f(2,bc%ibulk1) 
   da(3,:) = bc%invM2*f(3,bc%ibulk2)-bc%invM1*f(3,bc%ibulk1)
 
+  ! NOTE: In non-split nodes at fault edges M and f are assembled across the fault.
+  ! Hence, f1=f2, invM1=invM2=1/(M1+M2) instead of invMi=1/Mi, and da=0.
+ 
 end function get_weighted_jump
 
 !----------------------------------------------------------------------