[cig-commits] r6161 - in
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d: . includes
willic3 at geodynamics.org
willic3 at geodynamics.org
Fri Mar 2 11:06:28 PST 2007
Author: willic3
Date: 2007-03-02 11:06:27 -0800 (Fri, 02 Mar 2007)
New Revision: 6161
Added:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/fmin.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/includes/parmat_6.inc
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/invsymp.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lnsrch.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f
Log:
New and altered routines for new method of dealing with power-law
rheology. Things should compile, but testing/debugging is needed.
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am 2007-03-02 09:35:24 UTC (rev 6160)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am 2007-03-02 19:06:27 UTC (rev 6161)
@@ -54,6 +54,7 @@
elas_strs_mat_cmp_ss.F \
elastc.F \
fill.f \
+ fmin.f \
formdf_ss.f \
formes_ss.f \
formf_ss.f \
@@ -78,6 +79,7 @@
infellh.f \
infelqh.f \
invar.f \
+ invsymp.f \
iquate.f \
isort.f \
iterate.F \
@@ -89,6 +91,7 @@
lfit.f \
lflteq.f \
lnklst.f \
+ lnsrch.f \
load.f \
loadf.f \
loadx.f \
@@ -120,6 +123,7 @@
matmod_def.f \
meansh.f \
nchar.f \
+ newt.f \
nfind.f \
nnblnk.f \
open_append.F \
@@ -267,6 +271,7 @@
includes/nunits_dim.inc \
includes/nvisdat_def.inc \
includes/nvisdat_dim.inc \
+ includes/parmat_6.inc \
includes/prestr_matinit_ext.inc \
includes/prscal_dim.inc \
includes/rconsts.inc \
Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/fmin.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/fmin.f 2007-03-02 09:35:24 UTC (rev 6160)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/fmin.f 2007-03-02 19:06:27 UTC (rev 6161)
@@ -0,0 +1,72 @@
+c -*- Fortran -*-
+c
+c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c
+c PyLith by Charles A. Williams, Brad Aagaard, and Matt Knepley
+c
+c Copyright (c) 2004-2006 Rensselaer Polytechnic Institute
+c
+c Permission is hereby granted, free of charge, to any person obtaining
+c a copy of this software and associated documentation files (the
+c "Software"), to deal in the Software without restriction, including
+c without limitation the rights to use, copy, modify, merge, publish,
+c distribute, sublicense, and/or sell copies of the Software, and to
+c permit persons to whom the Software is furnished to do so, subject to
+c the following conditions:
+c
+c The above copyright notice and this permission notice shall be
+c included in all copies or substantial portions of the Software.
+c
+c THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+c EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+c MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+c NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+c LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+c OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+c WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+c
+c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c
+c
+ function fmin(x,n,fvec,rpar,nrpar,ipar,nipar,funcv,
+ & ierr,errstrng)
+c
+c... subroutine to compute half the dot product of a vector with
+c itself, given a function to compute the vector and associated
+c parameters.
+c
+ include "implicit.inc"
+c
+c... parameter definitions
+c
+ include "nconsts.inc"
+ include "rconsts.inc"
+c
+c... subroutine arguments
+c
+ integer n,nrpar,nipar,ierr
+ integer ipar(nipar)
+ double precision fmin,x(n),fvec(n),rpar(nrpar)
+ character errstrng*(*)
+ external funcv
+c
+c... external routines
+c
+ double precision ddot
+ external ddot
+c
+c... local variables
+c
+c
+ call funcv(n,x,fvec,rpar,nrpar,ipar,nipar,ierr,errstrng)
+ if(ierr.ne.izero) return
+ fmin=half*ddot(n,fvec,ione,fvec,ione)
+ return
+ end
+c
+c version
+c $Id: fmin.f,v 1.1 2004/04/14 21:18:30 willic3 Exp $
+c
+c Generated automatically by Fortran77Mill on Wed May 21 14:15:03 2003
+c
+c End of file
Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/includes/parmat_6.inc
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/includes/parmat_6.inc 2007-03-02 09:35:24 UTC (rev 6160)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/includes/parmat_6.inc 2007-03-02 19:06:27 UTC (rev 6161)
@@ -0,0 +1,61 @@
+c -*- Fortran -*-
+c
+c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c
+c PyLith by Charles A. Williams, Brad Aagaard, and Matt Knepley
+c
+c Copyright (c) 2004-2007 Rensselaer Polytechnic Institute
+c
+c Permission is hereby granted, free of charge, to any person obtaining
+c a copy of this software and associated documentation files (the
+c "Software"), to deal in the Software without restriction, including
+c without limitation the rights to use, copy, modify, merge, publish,
+c distribute, sublicense, and/or sell copies of the Software, and to
+c permit persons to whom the Software is furnished to do so, subject to
+c the following conditions:
+c
+c The above copyright notice and this permission notice shall be
+c included in all copies or substantial portions of the Software.
+c
+c THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+c EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+c MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+c NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+c LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+c OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+c WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+c
+c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c
+c... parmat_6.inc: Parameter definitions for material type 6.
+c
+c... The rpar array contains:
+c 1-6: Constant part of vector, which includes current total
+c strain, viscous strain from previous step, stress-related
+c quantities from previous time step, and prestress term.
+c 7-27: Inverted elastic material matrix.
+c 28: Integration parameter alfap.
+c 29: Time step size deltp.
+c 30: Viscosity coefficient emhu.
+c 31: Power-law exponent anpwr.
+c
+c... The ipar array contains:
+c 1-36: The iddmat index array.
+c
+ integer indconst,inddmati,indalfap,inddeltp,indemhu,indanpwr
+ integer indiddmat
+c
+ parameter(indconst=1,
+ & inddmati=7,
+ & indalfap=28,
+ & inddeltp=29,
+ & indemhu=30,
+ & indanpwr=31,
+ & indiddmat=1)
+c
+c version
+c $Id: parmat_6.inc,v 1.1 2004/04/14 21:18:30 willic3 Exp $
+c
+c Generated automatically by Fortran77Mill on Wed May 21 21:38:38 2003
+c
+c End of file
Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/invsymp.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/invsymp.f 2007-03-02 09:35:24 UTC (rev 6160)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/invsymp.f 2007-03-02 19:06:27 UTC (rev 6161)
@@ -0,0 +1,71 @@
+c -*- Fortran -*-
+c
+c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c
+c PyLith by Charles A. Williams, Brad Aagaard, and Matt Knepley
+c
+c Copyright (c) 2004-2007 Rensselaer Polytechnic Institute
+c
+c Permission is hereby granted, free of charge, to any person obtaining
+c a copy of this software and associated documentation files (the
+c "Software"), to deal in the Software without restriction, including
+c without limitation the rights to use, copy, modify, merge, publish,
+c distribute, sublicense, and/or sell copies of the Software, and to
+c permit persons to whom the Software is furnished to do so, subject to
+c the following conditions:
+c
+c The above copyright notice and this permission notice shall be
+c included in all copies or substantial portions of the Software.
+c
+c THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+c EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+c MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+c NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+c LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+c OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+c WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+c
+c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c
+c
+ subroutine invsymp(nsize,ndim,mat,ierr,errstrng)
+c
+c... routine to provide the inverse of a symmetric, positive definite
+c matrix stored in packed form.
+c This routine uses LAPACK routines.
+c
+ implicit none
+c
+c... parameter definitions
+c
+ include "nconsts.inc"
+c
+c... subroutine arguments
+c
+ integer nsize,ndim,ierr
+ double precision mat(nsize)
+ character errstrng*(*)
+c
+cdebug write(6,*) "Hello from invsymp_f!"
+c
+ call dpptrf('u',ndim,mat,ierr)
+ if(ierr.ne.izero) then
+ errstrng="invsymp (1)"
+ if(ierr.lt.izero) then
+ ierr=117
+ else
+ ierr=118
+ end if
+ return
+ end if
+ call dpptri('u',ndim,mat,ierr)
+ if(ierr.ne.izero) then
+ errstrng="invsymp (2)"
+ if(ierr.lt.izero) then
+ ierr=117
+ else
+ ierr=119
+ end if
+ end if
+ return
+ end
Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lnsrch.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lnsrch.f 2007-03-02 09:35:24 UTC (rev 6160)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lnsrch.f 2007-03-02 19:06:27 UTC (rev 6161)
@@ -0,0 +1,87 @@
+ SUBROUTINE lnsrch(n,xold,fold,g,p,x,f,stpmax,check,func,rpar,
+ & nrpar,ipar,nipar,funcv,ierr,errstrng)
+c
+c... routine to perform a line search along the direction given by p.
+c Adapted from Numerical Recipes.
+c
+ include "implicit.inc"
+c
+c... parameter definitions
+c
+ include "nconsts.inc"
+ include "rconsts.inc"
+ integer NP
+ double precision ALF,TOLX
+ PARAMETER (ALF=1.d-8,TOLX=1.d-12,NP=6)
+c
+c... subroutine arguments
+c
+ integer n,nrpar,nipar,ierr
+ integer ipar(nipar)
+ double precision xold(n),fold,g(n),p(n),x(n),f,stpmax,func
+ double precision rpar(nrpar)
+ external func,funcv
+ logical check
+ character errstrng*(*)
+CU USES func
+c
+c... external routines
+c
+ double precision dnrm2,ddot
+ external dnrm2,ddot
+c
+c... local variables
+c
+ INTEGER i
+ double precision a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2
+ double precision slope,sum,temp,test,tmplam,da
+ double precision fvec(NP)
+c
+ check=.false.
+ sum=dnrm2(n,p,ione)
+ if(sum.gt.stpmax)then
+ da=stpmax/sum
+ call dscal(n,da,p,ione)
+ endif
+ slope=ddot(n,g,ione,p,ione)
+ test=zero
+ do i=1,n
+ temp=abs(p(i))/max(abs(xold(i)),one)
+ if(temp.gt.test)test=temp
+ end do
+ alamin=TOLX/test
+ alam=one
+1 continue
+ do i=1,n
+ x(i)=xold(i)+alam*p(i)
+ end do
+ f=func(x,n,fvec,rpar,nrpar,ipar,nipar,funcv,ierr,errstrng)
+ if(alam.lt.alamin)then
+ call dcopy(n,xold,ione,x,ione)
+ check=.true.
+ return
+ else if(f.le.fold+ALF*alam*slope)then
+ return
+ else
+ if(alam.eq.one)then
+ tmplam=-slope/(two*(f-fold-slope))
+ else
+ rhs1=f-fold-alam*slope
+ rhs2=f2-fold2-alam2*slope
+ a=(rhs1/alam*alam-rhs2/(alam2*alam2))/(alam-alam2)
+ b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/(alam-alam2)
+ if(a.eq.zero)then
+ tmplam=-slope/(two*b)
+ else
+ disc=b*b-three*a*slope
+ tmplam=(-b+sqrt(disc))/(three*a)
+ endif
+ if(tmplam.gt.half*alam)tmplam=half*alam
+ endif
+ endif
+ alam2=alam
+ f2=f
+ fold2=fold
+ alam=max(tmplam,0.1d0*alam)
+ goto 1
+ END
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f 2007-03-02 09:35:24 UTC (rev 6160)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f 2007-03-02 19:06:27 UTC (rev 6161)
@@ -263,7 +263,6 @@
c
double precision e,pr,anpwr,emhu,rmu
double precision sdev(nstr),sinv1,steff
- double precision sigma,sigmap,sxx,syy,szz,sxy,syz,sxz
double precision cmat(nddmat)
c
c... included variable definitions
@@ -294,26 +293,12 @@
c Jacobian.
c
else
- call dpptrf('u',nstr,dmat,ierr)
- if(ierr.ne.izero) then
- errstrng="td_matinit_6(1)"
- if(ierr.lt.izero) then
- ierr=117
- else
- ierr=118
- end if
- end if
- call dpptri('u',nstr,dmat,ierr)
- if(ierr.ne.izero) then
- errstrng="td_matinit_6(2)"
- if(ierr.lt.izero) then
- ierr=117
- else
- ierr=119
- end if
- end if
- call fill(cmat,zero,nddmat)
c
+c... invert elastic matrix
+c
+ call invsymp(nddmat,nstr,dmat,ierr,errstrng)
+ if(ierr.ne.izero) return
+c
c... define material properties
c
e=prop(2)
@@ -322,70 +307,198 @@
emhu=prop(5)
rmu=half*e/(one+pr)
c
-c... get stress invariants and compute factors for Jacobian
+c... compute viscous Jacobian and add it to inverted elastic matrix
c
- sigma=(steff/emhu)**(anpwr-one)
- sigmap=(anpwr-one)*sigma
- sigma=half*alfap*deltp*sigma/emhu
- sigmap=fourth*alfap*deltp*sigmap/emhu
- sxx=sdev(1)/steff
- syy=sdev(2)/steff
- szz=sdev(3)/steff
- sxy=two*sdev(4)/steff
- syz=two*sdev(5)/steff
- sxz=two*sdev(6)/steff
+ call dbds_6(cmat,nddmat,iddmat,sdev,nstr,steff,anpwr,emhu,
+ & alfap,deltp)
+ call daxpy(nddmat,one,cmat,ione,dmat,ione)
c
-c... compute viscous Jacobian
+c... invert augmented matrix
c
- cmat(iddmat(1,1))=sigma*two*third+sigmap*sxx*sxx
- cmat(iddmat(1,2))=-sigma*third +sigmap*sxx*syy
- cmat(iddmat(1,3))=-sigma*third +sigmap*sxx*szz
- cmat(iddmat(1,4))= sigmap*sxx*sxy
- cmat(iddmat(1,5))= sigmap*sxx*syz
- cmat(iddmat(1,6))= sigmap*sxx*sxz
- cmat(iddmat(2,2))=sigma*two*third+sigmap*syy*syy
- cmat(iddmat(2,3))=-sigma*third +sigmap*syy*szz
- cmat(iddmat(2,4))= sigmap*syy*sxy
- cmat(iddmat(2,5))= sigmap*syy*syz
- cmat(iddmat(2,6))= sigmap*syy*sxz
- cmat(iddmat(3,3))=sigma*two*third+sigmap*szz*szz
- cmat(iddmat(3,4))= sigmap*szz*sxy
- cmat(iddmat(3,5))= sigmap*szz*syz
- cmat(iddmat(3,6))= sigmap*szz*sxz
- cmat(iddmat(4,4))=sigma*two +sigmap*sxy*sxy
- cmat(iddmat(4,5))= sigmap*sxy*syz
- cmat(iddmat(4,6))= sigmap*sxy*sxz
- cmat(iddmat(5,5))=sigma*two +sigmap*syz*syz
- cmat(iddmat(5,6))= sigmap*syz*sxz
- cmat(iddmat(6,6))=sigma*two +sigmap*sxz*sxz
+ call invsymp(nddmat,nstr,dmat,ierr,errstrng)
+ if(ierr.ne.izero) return
+ tmax=((emhu/steff)**(anpwr-one))*emhu/rmu
+ end if
+ return
+ end
c
-c... add this to inverted elastic matrix and invert again
c
- call daxpy(nddmat,one,cmat,ione,dmat,ione)
- call dpptrf('u',nstr,dmat,ierr)
- if(ierr.ne.izero) then
- errstrng="td_matinit_6(3)"
- if(ierr.lt.izero) then
- ierr=117
- else
- ierr=118
- end if
- end if
- call dpptri('u',nstr,dmat,ierr)
- if(ierr.ne.izero) then
- errstrng="td_matinit_6(4)"
- if(ierr.lt.izero) then
- ierr=117
- else
- ierr=119
- end if
- end if
- tmax=((emhu/steff)**(anpwr-one))*emhu/rmu
+ subroutine jaccmp_6(scur,n,fvec,NP,fjac,rpar,nrpar,ipar,nipar,
+ & ierr,errstrng)
+c
+c... subroutine to form the Jacobian used in finding the zero of the
+c stress function.
+c
+ include "implicit.inc"
+c
+c... parameter definitions
+c
+ include "ndimens.inc"
+ include "nconsts.inc"
+ include "rconsts.inc"
+ include "parmat_6.inc"
+c
+c... subroutine arguments
+c
+ integer n,NP,nrpar,nipar,ierr
+c*** note that the dimension below is a bit kludgy and does not allow
+c*** anything else to be put in the ipar array
+ integer ipar(nstr,nstr)
+ double precision scur(n),fvec(n),fjac(NP,NP),rpar(nrpar)
+ character errstrng*(*)
+c
+c... local variables
+c
+ integer i,j
+ctest integer iddmat(nstr,nstr)
+ctest equivalence(ipar(indiddmat),iddmat)
+ double precision dtmp(nddmat),cmat(nddmat)
+ double precision sdev(nstr),sinv1,steff
+c
+cdebug write(6,*) "Hello from jaccmp_6_f!"
+c
+c
+c... get stress invariants and a copy of inverted elasticity matrix
+c
+ call invar(sdev,sinv1,steff,scur)
+ call dcopy(nddmat,rpar(inddmati),ione,dtmp,ione)
+c
+c... if second deviatoric invariant is zero, use only elastic solution
+c
+ if(steff.eq.zero) then
+ call invsymp(nddmat,nstr,dtmp,ierr,errstrng)
+ if(ierr.ne.izero) return
+ do i=1,nstr
+ do j=1,nstr
+ fjac(i,j)=dtmp(ipar(i,j))
+ end do
+ end do
+ return
+c
+c... otherwise, invert elastic matrix and augment it by the viscous
+c Jacobian.
+c
+ else
+c
+c... compute viscous Jacobian and add it to inverted elastic matrix
+c
+ call dbds_6(cmat,nddmat,ipar,sdev,nstr,steff,
+ & rpar(indanpwr),rpar(indemhu),rpar(indalfap),rpar(inddeltp))
+ call daxpy(nddmat,one,cmat,ione,dtmp,ione)
+c
+c... invert augmented matrix and store results in fjac
+c
+ call invsymp(nddmat,nstr,dtmp,ierr,errstrng)
+ do i=1,nstr
+ do j=1,nstr
+ fjac(i,j)=dtmp(ipar(i,j))
+ end do
+ end do
end if
return
end
c
c
+ subroutine dbds_6(cmat,nddmat,iddmat,sdev,nstr,steff,anpwr,emhu,
+ & alfap,deltp)
+c
+c... subroutine to compute derivatives of viscous strain rate with
+c respect to stress.
+c
+ implicit none
+c
+c... parameter definitions
+c
+ include "nconsts.inc"
+ include "rconsts.inc"
+c
+c... subroutine arguments
+c
+ integer nddmat,nstr
+ integer iddmat(nstr,nstr)
+ double precision cmat(nddmat),sdev(nstr),steff,anpwr,emhu
+ double precision alfap,deltp
+c
+c... local variables
+c
+ double precision sigma,sigmap,sxx,syy,szz,sxy,syz,sxz
+c
+cdebug write(6,*) "Hello from dbds_6!"
+c
+ call fill(cmat,zero,nddmat)
+ sigma=(steff/emhu)**(anpwr-one)
+ sigmap=(anpwr-one)*sigma
+ sigma=half*alfap*deltp*sigma/emhu
+ sigmap=fourth*alfap*deltp*sigmap/emhu
+ sxx=sdev(1)/steff
+ syy=sdev(2)/steff
+ szz=sdev(3)/steff
+ sxy=two*sdev(4)/steff
+ syz=two*sdev(5)/steff
+ sxz=two*sdev(6)/steff
+c
+c... compute viscous Jacobian
+c
+ cmat(iddmat(1,1))=sigma*two*third+sigmap*sxx*sxx
+ cmat(iddmat(1,2))=-sigma*third +sigmap*sxx*syy
+ cmat(iddmat(1,3))=-sigma*third +sigmap*sxx*szz
+ cmat(iddmat(1,4))= sigmap*sxx*sxy
+ cmat(iddmat(1,5))= sigmap*sxx*syz
+ cmat(iddmat(1,6))= sigmap*sxx*sxz
+ cmat(iddmat(2,2))=sigma*two*third+sigmap*syy*syy
+ cmat(iddmat(2,3))=-sigma*third +sigmap*syy*szz
+ cmat(iddmat(2,4))= sigmap*syy*sxy
+ cmat(iddmat(2,5))= sigmap*syy*syz
+ cmat(iddmat(2,6))= sigmap*syy*sxz
+ cmat(iddmat(3,3))=sigma*two*third+sigmap*szz*szz
+ cmat(iddmat(3,4))= sigmap*szz*sxy
+ cmat(iddmat(3,5))= sigmap*szz*syz
+ cmat(iddmat(3,6))= sigmap*szz*sxz
+ cmat(iddmat(4,4))=sigma*two +sigmap*sxy*sxy
+ cmat(iddmat(4,5))= sigmap*sxy*syz
+ cmat(iddmat(4,6))= sigmap*sxy*sxz
+ cmat(iddmat(5,5))=sigma*two +sigmap*syz*syz
+ cmat(iddmat(5,6))= sigmap*syz*sxz
+ cmat(iddmat(6,6))=sigma*two +sigmap*sxz*sxz
+ return
+ end
+c
+c
+ subroutine betacmp_6(strs,beta,sdev,seff,sinv1,emhu,anpwr,nstr)
+c
+c... subroutine to compute viscous strain rate vector beta
+c Routine also returns stress invariants
+c
+ implicit none
+c
+c... parameter definitions
+c
+ include "nconsts.inc"
+ include "rconsts.inc"
+c
+c... subroutine arguments
+c
+ integer nstr
+ double precision strs(nstr),beta(nstr),sdev(nstr),seff,sinv1
+ double precision emhu,anpwr
+c
+c... local variables
+c
+ double precision rm
+c
+c... compute stress invariants and constants for loop
+c
+ call fill(beta,zero,nstr)
+ call invar(sdev,sinv1,seff,strs)
+ rm=half*(seff/emhu)**(anpwr-one)/emhu
+c
+c... compute strain rate components
+c
+ call daxpy(nstr,rm,sdev,ione,beta,ione)
+ return
+ end
+c
+c
subroutine td_strs_6(state,dstate,state0,ee,scur,dmat,prop,
& rtimdat,rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,
& matchg,ierr,errstrng)
@@ -402,6 +515,9 @@
include "ndimens.inc"
include "nconsts.inc"
include "rconsts.inc"
+ include "parmat_6.inc"
+ integer nrpar,nipar
+ parameter(nrpar=31,nipar=36)
c
c... subroutine arguments
c
@@ -420,21 +536,22 @@
c
c... external functions
c
- double precision sprod
- external sprod
+ double precision sprod,dasum
+ external sprod,funcv_6,jaccmp_6,dasum
c
c... local constants
c
+ double precision sub
+ data sub/-1.0d0/
c
c... local variables
c
integer i
double precision e,pr,anpwr,emhu,rmu
- double precision sinv1t,sefft,sinv1tdt,sefftdt,beta,betat,betatdt
- double precision rmt,rmtdt
- double precision sdevt(nstr),sdevtdt(nstr),eetdt(nstr)
- double precision dtmp(nddmat)
-
+ double precision test0,rm,rmt,rmtdt
+ double precision sdev(nstr),betat(nstr),betatdt(nstr),sinv1,seff
+ double precision rpar(nrpar)
+ logical check
c
c... included variable definitions
c
@@ -445,6 +562,7 @@
cdebug write(6,*) "Hello from td_strs_6!"
c
ierr=0
+ test0=dasum(nstate0,state0,ione)
c
c... define material properties
c
@@ -454,48 +572,124 @@
emhu=prop(5)
rmu=half*e/(one+pr)
tmax=big
+ rmt=deltp*(one-alfap)
+ rm=sub*rmt
+ rmtdt=deltp*alfap
c
-c... for first iteration, beta is computed using only the stresses at
-c time t. Otherwise, the stress estimates at time t + dt are also
-c used.
+c... for first iteration, use stresses from previous step as initial
+c guess, otherwise use current estimate.
c
- call invar(sdevt,sinv1t,sefft,state)
- if(iter.eq.ione) then
- do i=1,nstr
- beta=half*(sefft/emhu)**(anpwr-one)*sdevt(i)/emhu
- dstate(i+12)=state(i+12)+deltp*beta
- eetdt(i)=ee(i)-dstate(i+12)
- end do
- else
- rmt=one-alfap
- rmtdt=alfap
- call invar(sdevtdt,sinv1tdt,sefftdt,dstate)
- do i=1,nstr
- betat=half*(sefft/emhu)**(anpwr-one)*sdevt(i)/emhu
- betatdt=half*(sefftdt/emhu)**(anpwr-one)*sdevtdt(i)/emhu
- dstate(i+12)=state(i+12)+deltp*(rmt*betat+rmtdt*betatdt)
- eetdt(i)=ee(i)-dstate(i+12)
- end do
- end if
+ if(iter.eq.ione) call dcopy(nstr,state,ione,dstate,ione)
c
-c... compute current stress estimates and copy state variables into the
-c appropriate spot
+c... compute constant part of iterative solution equation.
c
- call elas_mat_6(dtmp,prop,iddmat,nprop,ierr,errstrng)
+c... current total strain
+ call dcopy(nstr,ee,ione,rpar(indconst),ione)
+c... inverted elasticity matrix times initial stresses
+ call elas_mat_6(rpar(inddmati),prop,iddmat,nprop,ierr,errstrng)
+ call invsymp(nddmat,nstr,rpar(inddmati),ierr,errstrng)
+ if(ierr.ne.izero) return
+ if(test0.ne.zero) call dspmv("u",nstr,one,rpar(inddmati),state0,
+ & ione,one,rpar(indconst),ione)
+c... viscous strain from previous time step
+ call daxpy(nstr,sub,state(13),ione,rpar(indconst),ione)
+c... contribution to viscous strain from stresses at beginning of step
+ call betacmp_6(state,betat,sdev,seff,sinv1,emhu,anpwr,nstr)
+ call daxpy(nstr,rm,betat,ione,rpar(indconst),ione)
+c
+c... finish defining parameters for stress computation then call Newton
+c routine to compute stresses.
+c
+ rpar(indalfap)=alfap
+ rpar(inddeltp)=deltp
+ rpar(indemhu)=emhu
+ rpar(indanpwr)=anpwr
+ call newt(dstate,nstr,rpar,nrpar,iddmat,nipar,funcv_6,jaccmp_6,
+ & check,ierr,errstrng)
+ if(ierr.ne.izero) return
+c
+c... update state variables for current stress estimates.
+c
+ call betacmp_6(dstate,betatdt,sdev,seff,sinv1,emhu,anpwr,nstr)
call dcopy(nstr,ee,ione,dstate(7),ione)
- call dcopy(nstr,state0,ione,dstate,ione)
- call dspmv("u",nstr,one,dtmp,eetdt,ione,one,dstate,ione)
+ call dcopy(nstr,state(13),ione,dstate(13),ione)
+ call daxpy(nstr,rmt,betat,ione,dstate(13),ione)
+ call daxpy(nstr,rmtdt,betatdt,ione,dstate(13),ione)
call dcopy(nstr,dstate,ione,scur,ione)
c
c... compute current Maxwell time
c
- call invar(sdevtdt,sinv1tdt,sefftdt,dstate)
- if(sefftdt.ne.zero) tmax=((emhu/sefftdt)**(anpwr-one))*emhu/rmu
+ call invar(sdev,sinv1,seff,dstate)
+ if(seff.ne.zero) tmax=((emhu/seff)**(anpwr-one))*emhu/rmu
c
return
end
c
c
+ subroutine funcv_6(n,scur,r,rpar,nrpar,ipar,nipar,
+ & ierr,errstrng)
+c
+c... routine to compute the vector of functions(r) evaluated at the
+c given stress state (contained in scur).
+c
+c... The rpar array contains:
+c 1-6: Constant part of vector, which includes current total
+c strain, viscous strain from previous step, stress-related
+c quantities from prefious time step, and prestress term.
+c 7-27: Inverted elastic material matrix.
+c 28: Integration parameter alfap.
+c 29: Time step size deltp.
+c 30: Viscosity coefficient emhu.
+c 31: Power-law exponent anpwr.
+c
+c... The ipar array contains:
+c 1-36: The iddmat index array
+c
+ implicit none
+c
+c... parameter definitions
+c
+ include "ndimens.inc"
+ include "nconsts.inc"
+ include "rconsts.inc"
+ include "parmat_6.inc"
+c
+c... subroutine arguments
+c
+ integer n,nrpar,nipar,ierr
+ integer ipar(nipar)
+ double precision scur(n),r(n),rpar(nrpar)
+ character errstrng*(*)
+c
+c... local data
+c
+ double precision sub
+ data sub/-1.0d0/
+c
+c... local variables
+c
+ double precision betatdt(nstr),rmtdt
+ double precision sdev(nstr),seff,sinv1
+c
+c... viscous strain rate multiplier
+c
+ rmtdt=-rpar(indalfap)*rpar(inddeltp)
+c
+c... compute current viscous strain rate
+c
+ call betacmp_6(scur,betatdt,sdev,seff,sinv1,
+ & rpar(indemhu),rpar(indanpwr),nstr)
+c
+c... compute contributions to r and accumulate
+c
+ call dcopy(nstr,rpar(indconst),ione,r,ione)
+ call daxpy(nstr,rmtdt,betatdt,ione,r,ione)
+ call dspmv("u",nstr,sub,rpar(inddmati),scur,ione,one,r,ione)
+c
+ return
+ end
+c
+c
subroutine td_strs_mat_6(state,dstate,state0,ee,scur,dmat,prop,
& rtimdat,rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,
& matchg,ierr,errstrng)
Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f 2007-03-02 09:35:24 UTC (rev 6160)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f 2007-03-02 19:06:27 UTC (rev 6161)
@@ -0,0 +1,105 @@
+ SUBROUTINE newt(x,n,rpar,nrpar,ipar,nipar,funcv,jaccmp,check,
+ & ierr,errstrng)
+c
+c... routine to find zero of a function using Newton's method with
+c line search.
+c Adapted from Numerical Recipes.
+c
+ include "implicit.inc"
+c
+c... parameter definitions
+c
+ include "nconsts.inc"
+ include "rconsts.inc"
+ integer NP,MAXITS
+ double precision TOLF,TOLMIN,TOLX,STPMX
+ parameter(NP=6,MAXITS=100,TOLF=1.0d-8,TOLMIN=1.0d-10,TOLX=1.0d-12,
+ & STPMX=100.0d0)
+c
+c... subroutine arguments
+c
+ integer n,nrpar,nipar,ierr
+ integer ipar(nipar)
+ double precision x(n),rpar(nrpar)
+ logical check
+ character errstrng*(*)
+c
+c... external routines
+c
+ double precision ddot,dnrm2
+ external funcv,jaccmp,ddot,dnrm2
+c
+c... local variables
+c
+ double precision fvec(NP)
+CU USES fdjac,fmin,lnsrch,lubksb,ludcmp
+ INTEGER i,its,j,indx(NP)
+ double precision d,den,f,fold,stpmax,sum,temp,test,fjac(NP,NP)
+ double precision g(NP),p(NP),xold(NP),fmin
+ external fmin
+c
+ f=fmin(x,n,fvec,rpar,nrpar,ipar,nipar,funcv,ierr,errstrng)
+ if(ierr.ne.izero) return
+ test=zero
+ do i=1,n
+ if(abs(fvec(i)).gt.test)test=abs(fvec(i))
+ end do
+ if(test.lt.0.01d0*TOLF)return
+ sum=dnrm2(n,x,ione)
+ stpmax=STPMX*max(sum,dble(n))
+ do its=1,MAXITS
+ call jaccmp(x,n,fvec,NP,fjac,rpar,nrpar,ipar,nipar,
+ & ierr,errstrng)
+c**** I should replace this with a BLAS call
+ do i=1,n
+ sum=zero
+ do j=1,n
+ sum=sum+fjac(j,i)*fvec(j)
+ end do
+ g(i)=sum
+ end do
+ call dcopy(n,x,ione,xold,ione)
+ fold=f
+ do i=1,n
+ p(i)=-fvec(i)
+ end do
+c**** these can be replaced by LAPACK calls -- eventually, it will
+c**** probably make more sense to use symmetric packed format for
+c**** all matrices.
+ call ludcmp(fjac,n,NP,indx,d)
+ call lubksb(fjac,n,NP,indx,p)
+ call lnsrch(n,xold,fold,g,p,x,f,stpmax,check,fmin,rpar,nrpar,
+ & ipar,nipar,funcv,ierr,errstrng)
+ if(ierr.ne.izero) return
+ test=zero
+ do i=1,n
+ if(abs(fvec(i)).gt.test)test=abs(fvec(i))
+ end do
+ if(test.lt.TOLF)then
+ check=.false.
+ return
+ endif
+ if(check)then
+ test=zero
+ den=max(f,half*dble(n))
+ do i=1,n
+ temp=abs(g(i))*max(abs(x(i)),one)/den
+ if(temp.gt.test)test=temp
+ end do
+ if(test.lt.TOLMIN)then
+ check=.true.
+ else
+ check=.false.
+ endif
+ return
+ endif
+ test=zero
+ do i=1,n
+ temp=(abs(x(i)-xold(i)))/max(abs(x(i)),one)
+ if(temp.gt.test)test=temp
+ end do
+ if(test.lt.TOLX)return
+ end do
+c**** replace this with an error code
+ pause 'MAXITS exceeded in newt'
+ END
More information about the cig-commits
mailing list