Removed passing of call_cnt for debugging, extra debugging output, added...

Removed passing of call_cnt for debugging, extra debugging output, added option of using 2D strain rate invariants and fixed Mohr-Coulomb plasticity formatting to be similar to the other options

src/vrm.f90

@@ -18,7 +18,7 @@
 
 subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure,             &
                plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb,   &
-               straintot,call_cnt)
+               straintot)
 
 use constants
 use definitions
@@ -42,7 +42,6 @@ double precision :: J2d,J3d,viscosity,pressure,compressibility,straintot
 character (len=8) plasticity_type
 double precision plasticity_parameters(9)
 logical ::  is_plastic,flag_vrm_pb
-integer call_cnt
 
 !------------------------------------------------------------------------------|
 !(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
@@ -57,14 +56,6 @@ double precision :: strain_soft_yield
 
 flag_vrm_pb = .false.
 
-!write (*,*) 'vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv'
-!write (*,*) 'init_e2d: ',params%init_e2d
-!write (*,*) 'real e2d: ',sqrt(J2d)
-!write (*,*) 'viscosity: ',viscosity
-!write (*,*) 'is_plastic: ',is_plastic
-!write (*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
-!write (*,*) ''
-
 if (params%init_e2d) then
   e2d=params%e2d0
 else
@@ -111,17 +102,6 @@ select case (trim(plasticity_type))
     k     = (6.d0*c*cos(phi))/(sqrt(3.d0)*(3.d0-sin(phi)))
     yield = max(c/compressibility,k-alpha*pressure)
     fail  = 2.d0*viscosity*e2d-yield
-    !if (call_cnt.le.10) then
-    !  write(*,*) '******************************************'
-    !  write(*,*) 'phi: ',phi
-    !  write(*,*) 'c: ',c
-    !  write(*,*) 'alpha: ',alpha
-    !  write(*,*) 'k: ',k
-    !  write(*,*) 'pressure: ',pressure
-    !  write(*,*) 'yield: ',yield
-    !  write(*,*) 'fail: ',fail
-    !  write(*,*) '******************************************'
-    !endif
     if (fail.gt.0.d0) then
       viscosity=0.5d0*yield/e2d
       is_plastic=.true.
@@ -287,16 +267,11 @@ select case (trim(plasticity_type))
     c = plasticity_parameters(2)
     sin_phi = sin(phi)
     cos_phi = cos(phi)
-    zeta      = cos_theta-sin_phi*sin_theta/sqrt3
-    ! modified by Jean
-    !fail      = -pressure*sin_phi+e2d*zeta-c*cos_phi
-    fail      = pressure*sin_phi+2.d0*viscosity*e2d*zeta-c*cos_phi
+    zeta    = cos_theta-sin_phi*sin_theta/sqrt3
+    yield = max(plasticity_parameters(2)/compressibility,(-pressure*sin_phi+c*cos_phi)/zeta)
+    fail  = 2.d0*viscosity*e2d-yield
     if (fail.gt.0.d0) then
-      viscosity=0.5d0/e2d*(c*cos_phi-pressure*sin_phi)/zeta
-      !if (viscosity<0.d0) then
-      !  viscosity=0.d0
-      !  !print *,viscosity,c*cos_phi,pressure*sin_phi
-      !end if
+      viscosity=0.5d0*yield/e2d
       is_plastic=.true.
     else
       is_plastic=.false.
@@ -325,14 +300,6 @@ select case (trim(plasticity_type))
 
 end select
 
-!write (*,*) 'vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv'
-!write (*,*) 'init_e2d: ',params%init_e2d
-!write (*,*) 'e2d out: ',sqrt(J2d)
-!write (*,*) 'viscosity out: ',viscosity
-!write (*,*) 'is_plastic out: ',is_plastic
-!write (*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
-!write (*,*) ''
-
 return
 end
 !-------------------------------------------------------------------------------