[cig-commits] commit: add a control of the rake of afterslip, export slip model to gmt and tune time steps.
Mercurial
hg at geodynamics.org
Tue Sep 20 12:13:00 PDT 2011
changeset: 5:21d47e0a9bf4
user: Sylvain Barbot <sylbar.vainbot at gmail.com>
date: Mon Apr 11 08:30:06 2011 -0700
files: elastic3d.f90 export.f90 friction3d.f90 include.f90 relax.f90
description:
add a control of the rake of afterslip, export slip model to gmt and tune time steps.
diff -r 2329ee35e5fb -r 21d47e0a9bf4 elastic3d.f90
--- a/elastic3d.f90 Thu Jan 06 18:12:33 2011 -0800
+++ b/elastic3d.f90 Mon Apr 11 08:30:06 2011 -0700
@@ -37,7 +37,7 @@ MODULE elastic3d
TYPE PLANE_STRUCT
SEQUENCE
- REAL*8 :: x,y,z,width,length,strike,dip
+ REAL*8 :: x,y,z,width,length,strike,dip,rake
END TYPE PLANE_STRUCT
TYPE LAYER_STRUCT
@@ -73,7 +73,7 @@ MODULE elastic3d
TYPE EVENT_STRUC
REAL*8 :: time
- INTEGER*4 :: ns,nt,nm,nl
+ INTEGER*4 :: i,ns,nt,nm,nl
TYPE(SOURCE_STRUCT), DIMENSION(:), ALLOCATABLE :: s,sc,ts,tsc,m,mc,l,lc
END TYPE EVENT_STRUC
diff -r 2329ee35e5fb -r 21d47e0a9bf4 export.f90
--- a/export.f90 Thu Jan 06 18:12:33 2011 -0800
+++ b/export.f90 Mon Apr 11 08:30:06 2011 -0700
@@ -733,7 +733,7 @@ END SUBROUTINE exporteigenstrain
DO k=1,np
CALL monitorfriction(n(k)%x,n(k)%y,n(k)%z, &
- n(k)%width,n(k)%length,n(k)%strike,n(k)%dip,beta, &
+ n(k)%width,n(k)%length,n(k)%strike,n(k)%dip,n(k)%rake,beta, &
sx1,sx2,sx3,dx1,dx2,dx3,sig,structure,slippatch)
ns1=SIZE(slippatch,1)
@@ -1022,6 +1022,83 @@ END SUBROUTINE exportcreep
CLOSE(15)
END SUBROUTINE exportvtk_grid
+
+ !------------------------------------------------------------------
+ ! subroutine ExportXY_RFaults
+ ! creates a .xy file (in the GMT closed-polygon format) containing
+ ! the rectangular faults. Each fault segemnt is described by a
+ ! closed polygon (rectangle) associated with a slip amplitude.
+ ! use pxzy with the -Cpalette.cpt -L -M options to color rectangles
+ ! by slip.
+ !
+ ! sylvain barbot 03/05/11 - original form
+ !------------------------------------------------------------------
+ SUBROUTINE exportxy_rfaults(e,rffilename)
+ TYPE(EVENT_STRUC), INTENT(IN) :: e
+ CHARACTER(80), INTENT(IN) :: rffilename
+
+ INTEGER :: iostatus,k
+ CHARACTER :: q
+
+ REAL*8 :: strike,dip,x1,x2,x3,cstrike,sstrike,cdip,sdip,L,W,slip
+
+ REAL*8, DIMENSION(3) :: s,d
+
+ ! double-quote character
+ q=char(34)
+
+ OPEN (UNIT=15,FILE=rffilename,IOSTAT=iostatus,FORM="FORMATTED")
+ IF (iostatus>0) THEN
+ WRITE_DEBUG_INFO
+ PRINT '(a)', rffilename
+ STOP "could not open file for export"
+ END IF
+
+ DO k=1,e%ns
+
+ ! fault slip
+ slip=e%s(k)%slip
+
+ ! fault orientation
+ strike=e%s(k)%strike
+ dip=e%s(k)%dip
+
+ ! fault center position
+ x1=e%s(k)%x
+ x2=e%s(k)%y
+ x3=e%s(k)%z
+
+ ! fault dimension
+ W=e%s(k)%width
+ L=e%s(k)%length
+
+ cstrike=cos(strike)
+ sstrike=sin(strike)
+ cdip=cos(dip)
+ sdip=sin(dip)
+
+ ! strike-slip unit direction
+ s(1)=sstrike
+ s(2)=cstrike
+ s(3)=0._8
+
+ ! dip-slip unit direction
+ d(1)=+cstrike*sdip
+ d(2)=-sstrike*sdip
+ d(3)=+cdip
+
+ ! fault edge coordinates
+ WRITE (15,'("> -Z",3ES11.2)') ABS(slip)
+ WRITE (15,'(3ES11.2)') x1-d(1)*W/2-s(1)*L/2, x2-d(2)*W/2-s(2)*L/2
+ WRITE (15,'(3ES11.2)') x1-d(1)*W/2+s(1)*L/2, x2-d(2)*W/2+s(2)*L/2
+ WRITE (15,'(3ES11.2)') x1+d(1)*W/2+s(1)*L/2, x2+d(2)*W/2+s(2)*L/2
+ WRITE (15,'(3ES11.2)') x1+d(1)*W/2-s(1)*L/2, x2+d(2)*W/2-s(2)*L/2
+
+ END DO
+
+ CLOSE(15)
+
+ END SUBROUTINE exportxy_rfaults
!------------------------------------------------------------------
! subroutine ExportVTK_RFaults
diff -r 2329ee35e5fb -r 21d47e0a9bf4 friction3d.f90
--- a/friction3d.f90 Thu Jan 06 18:12:33 2011 -0800
+++ b/friction3d.f90 Mon Apr 11 08:30:06 2011 -0700
@@ -33,32 +33,35 @@ CONTAINS
!-----------------------------------------------------------------
! subroutine FrictionPlaneExpEigenStress
+ !
+ ! *** this function is deprecated ***
+ !
! compute the eigen-stress (forcing moment) to be relaxed by
! rate-dependent inelastic deformation in the case of a frictional
! surface:
!
- ! sigma^i = C:F:sigma
+ ! sigma^i = C:F:sigma
!
! where C is the elastic moduli tensor, F is the heterogeneous
! fluidity moduli tensor and sigma is the instantaneous stress
! tensor. for a frictional surface, the eigenstrain-rate is given
! by
!
- ! epsilon^i^dot = F:sigma = gamma^dot R
+ ! epsilon^i^dot = F:sigma = gamma^dot R
!
! where gamma^dot is the slip rate (a scalar) and R is the
! deviatoric, symmetric, and unitary, tensor:
!
- ! R_ij = 1/2 ( t_i n_j + t_j n_i )
+ ! R_ij = 1/2 ( t_i n_j + t_j n_i ) / sqrt( t_j t_j )
!
! where the shear traction t_i is the projection of the traction
! vector on the plane surface. the strain amplitude is given by
!
- ! gamma^dot = 2 vo sinh( taus / (t_c )
+ ! gamma^dot = vo sinh( taus / (t_c )
!
! where taus is the effective shear on the fault plane,
!
- ! taus = tau + mu*sigma
+ ! taus = tau + mu*sigma
!
! where tau is the shear and sigma the normal stress. tau and sigma
! assumed to be the co-seismic change only, not the absolute
@@ -66,12 +69,13 @@ CONTAINS
! stress, corresponds to (a-b)*sigma in the framework of rate-and-
! state friction. the effective viscosity eta* and the fluidity
!
- ! eta* = tau / gamma^dot
+ ! eta* = tau / gamma^dot
! fluidity = 1 / eta*
!
- ! are used to compute the optimal time-step.
+ ! are used to compute the optimal time-step.
!
! sylvain barbot (07/24/07) - original form
+ ! (07/24/07) - deprecated (see frictioneigenstress)
!-----------------------------------------------------------------
SUBROUTINE frictionplaneeigenstress(sig,mu,structure, &
n1,n2,n3,sx1,sx2,sx3,dx1,dx2,dx3,moment,maxwelltime,gamma,dt)
@@ -187,21 +191,21 @@ CONTAINS
! tensor. for a frictional surface, the eigenstrain-rate is given
! by
!
- ! epsilon^i^dot = F:sigma = gamma^dot R
+ ! epsilon^i^dot = F:sigma = gamma^dot R
!
! where gamma^dot is the slip rate (a scalar) and R is the
! deviatoric, symmetric, and unitary, tensor:
!
- ! R_ij = 1/2 ( t_i n_j + t_j n_i )
+ ! R_ij = 1/2 ( t_i n_j + t_j n_i ) / sqrt( t_j t_j )
!
! where the shear traction t_i is the projection of the traction
! vector on the plane surface. the strain amplitude is given by
!
- ! gamma^dot = 2 vo sinh( taus / (t_c )
+ ! gamma^dot = H( t_j r_j ) 2 vo sinh( taus / (t_c )
!
! where taus is the effective shear on the fault plane,
!
- ! taus = tau + mu*sigma
+ ! taus = tau + mu*sigma
!
! where tau is the shear and sigma the normal stress. tau and sigma
! assumed to be the co-seismic change only, not the absolute
@@ -209,17 +213,22 @@ CONTAINS
! stress, corresponds to (a-b)*sigma in the framework of rate-and-
! state friction. the effective viscosity eta* and the fluidity
!
- ! eta* = tau / gamma^dot
- ! fluidity = 1 / eta*
+ ! eta* = tau / gamma^dot
+ ! fluidity = 1 / eta*
!
- ! are used to compute the optimal time-step.
+ ! are used to compute the optimal time-step. H( x ) is the
+ ! Heaviside function and r_i is the rake vector. I impose
+ ! gamma^dot to be zero is t_j r_j < 0. This constraint is
+ ! enforced to ensure that no back slip occurs on faults.
!
! sylvain barbot (07/24/07) - original form
+ ! (02/28/11) - add constraints on the direction of
+ ! afterslip
!-----------------------------------------------------------------
- SUBROUTINE frictioneigenstress(x,y,z,L,W,strike,dip,beta, &
+ SUBROUTINE frictioneigenstress(x,y,z,L,W,strike,dip,rake,beta, &
sig,mu,structure,sx1,sx2,sx3,dx1,dx2,dx3,moment,maxwelltime,vel)
INTEGER, INTENT(IN) :: sx1,sx2,sx3
- REAL*8, INTENT(IN) :: mu,dx1,dx2,dx3,x,y,z,L,W,strike,dip,beta
+ REAL*8, INTENT(IN) :: mu,dx1,dx2,dx3,x,y,z,L,W,strike,dip,rake,beta
TYPE(LAYER_STRUCT), DIMENSION(:), INTENT(IN) :: structure
TYPE(TENSOR), INTENT(IN), DIMENSION(sx1,sx2,sx3) :: sig
TYPE(TENSOR), INTENT(INOUT), DIMENSION(sx1,sx2,sx3) :: moment
@@ -232,11 +241,11 @@ CONTAINS
INTEGER :: i1,i2,i3
TYPE(TENSOR) :: s
- REAL*8, DIMENSION(3) :: t,ts,n
+ REAL*8, DIMENSION(3) :: t,ts,n,r
REAL*8 :: vo,tauc,taun,taus,gammadot,impulse, &
friction,tau,scaling,cohesion
REAL*8 :: x1,x2,x3,x1s,x2s,x3s,x1i,x3i, &
- cstrike,sstrike,cdip,sdip,x2r,&
+ cstrike,sstrike,cdip,sdip,cr,sr,x2r,&
temp1,temp2,temp3,sourc,image,xr,yr,zr,Wp,Lp,dum
REAL*4 :: tm
@@ -253,6 +262,8 @@ CONTAINS
sstrike=sin(strike)
cdip=cos(dip)
sdip=sin(dip)
+ cr=cos(rake)
+ sr=sin(rake)
! effective tapered dimensions
Wp=W*(1._8+2._8*beta)/2._8
@@ -269,6 +280,11 @@ CONTAINS
n(2)=-cdip*sstrike
n(3)=-sdip
+ ! rake vector component
+ r(1)=sstrike*cr+cstrike*sdip*sr
+ r(2)=cstrike*cr-sstrike*sdip*sr
+ r(3)=+cdip*sr
+
DO i3=1,sx3
x3=DBLE(i3-1)*dx3
IF ((abs(x3-z).gt.Lp) .and. (abs(x3+z).gt.Lp)) CYCLE
@@ -292,8 +308,8 @@ CONTAINS
x3s= sdip*x2r+cdip*x3
x3i=-sdip*x2r+cdip*x3
- !integrate at depth and along strike with raised cosine taper
- !and shift sources to x,y,z coordinate
+ ! integrate at depth and along strike with raised cosine taper
+ ! and shift sources to x,y,z coordinate
temp1=gauss(x1s-xr,dx1)
temp2=omega((x2s-yr)/W,beta)
temp3=omega((x3s-zr)/L,beta)
@@ -318,6 +334,9 @@ CONTAINS
! effective shear stress on fault plane
tau=taus+friction*taun-cohesion
+
+ ! rake direction test
+ IF (SUM(ts*r).LT.0.d0) CYCLE
! warning for wrong input
IF ((tau/tauc) .gt. 20) THEN
@@ -381,26 +400,28 @@ CONTAINS
!
! sylvain barbot (10-16-07) - original form
!---------------------------------------------------------------------
- SUBROUTINE monitorfriction(x,y,z,L,W,strike,dip,beta, &
+ SUBROUTINE monitorfriction(x,y,z,L,W,strike,dip,rake,beta, &
sx1,sx2,sx3,dx1,dx2,dx3,sig,structure,patch)
INTEGER, INTENT(IN) :: sx1,sx2,sx3
- REAL*8, INTENT(IN) :: x,y,z,L,W,strike,dip,beta,dx1,dx2,dx3
+ REAL*8, INTENT(IN) :: x,y,z,L,W,strike,rake,dip,beta,dx1,dx2,dx3
TYPE(TENSOR), DIMENSION(sx1,sx2,sx3), INTENT(IN) :: sig
TYPE(SLIPPATCH_STRUCT), ALLOCATABLE, DIMENSION(:,:), INTENT(OUT) :: patch
TYPE(LAYER_STRUCT), DIMENSION(:), INTENT(IN) :: structure
INTEGER :: i1,i2,i3,px2,px3,j2,j3,status
- REAL*8 :: cstrike,sstrike,cdip,sdip,slip,ss,ds
+ REAL*8 :: cstrike,sstrike,cdip,sdip,slip,ss,ds,cr,sr
REAL*8 :: vo,tauc,taun,taus, &
friction,tau,cohesion
REAL*8 :: x1,x2,x3,xr,yr,zr,Wp,Lp
TYPE(TENSOR) :: s
- REAL*8, DIMENSION(3) :: t,ts,n,sv,dv
+ REAL*8, DIMENSION(3) :: t,ts,n,sv,dv,r
cstrike=cos(strike)
sstrike=sin(strike)
cdip=cos(dip)
sdip=sin(dip)
+ cr=cos(rake)
+ sr=sin(rake)
! strike direction vector
sv=(/ sstrike, cstrike, 0._8 /)
@@ -425,6 +446,11 @@ CONTAINS
n(2)=-cdip*sstrike
n(3)=-sdip
+ ! rake vector component
+ r(1)=sstrike*cr+cstrike*sdip*sr
+ r(2)=cstrike*cr-sstrike*sdip*sr
+ r(3)=+cdip*sr
+
! loop in the dip direction
DO j3=1,px3+1
! loop in the strike direction
@@ -469,13 +495,15 @@ CONTAINS
! effective shear stress on fault plane
tau=taus+friction*taun-cohesion
- ! yield surface test
- IF ((0._8 .GE. taus) .OR. (tau .LE. 0._8)) THEN
+ ! shear traction direction
+ ts=ts/taus
+
+ ! yield surface rake direction tests
+ IF ((0._8 .GE. taus) .OR. &
+ (tau .LE. 0._8) .OR. &
+ (SUM(ts*r).LT.0.d0)) THEN
ss=0;ds=0;slip=0;
ELSE
- ! shear traction direction
- ts=ts/taus
-
! creep rate
slip=vo*2._8*sinh(tau/tauc)
diff -r 2329ee35e5fb -r 21d47e0a9bf4 include.f90
--- a/include.f90 Thu Jan 06 18:12:33 2011 -0800
+++ b/include.f90 Mon Apr 11 08:30:06 2011 -0700
@@ -21,11 +21,11 @@
! export amplitude of scalar fields
! along a plane in GRD binary format
-!#define GRD_EXPORTEIGENSTRAIN 1
+#define GRD_EXPORTEIGENSTRAIN 1
! export creep velocity along a frictional
! plane in GRD binary format
-!#define GRD_EXPORTCREEP 1
+#define GRD_EXPORTCREEP 1
! export data to the TXT format
!#define TXT 1
@@ -35,7 +35,7 @@
! export amplitude of scalar fields along
! an observation plane in text format
-!#define TXT_EXPORTEIGENSTRAIN 1
+#define TXT_EXPORTEIGENSTRAIN 1
! export creep velocity along a frictional
! plane in text format
diff -r 2329ee35e5fb -r 21d47e0a9bf4 relax.f90
--- a/relax.f90 Thu Jan 06 18:12:33 2011 -0800
+++ b/relax.f90 Mon Apr 11 08:30:06 2011 -0700
@@ -72,7 +72,7 @@ PROGRAM relax
! N (x1)
! /
! /| Strike
- ! Pos:-> @------------------------ (x2)
+ ! x1,x2,x3 ->@------------------------ (x2)
! |\ p . \ W
! :-\ i . \ i
! | \ l . \ d
@@ -86,6 +86,11 @@ PROGRAM relax
! Dislocations are converted to double-couple equivalent body-force
! analytically. Solution displacement is obtained by application of
! the Greens functions in the Fourier domain.
+ !
+ ! For friction faults where slip rates are evaluated from stress and
+ ! a constitutive law, the rake corresponds to the orientation of slip.
+ ! That is, if r_i is the rake vector and v_i is the instantaneous
+ ! velocity vector, then r_j v_j >= 0.
!
! OUTPUT:
! The vector-valued deformation is computed everywhere in a cartesian
@@ -148,6 +153,9 @@ PROGRAM relax
! and output components of stress tensor
! (07-19-10) - includes surface tractions initial condition
! output geometry in VTK format for Paraview
+ ! (02-28-11) - add constraints on the broad direction of
+ ! afterslip, export faults to GMT xy format
+ ! and allow scaling of computed time steps.
!-----------------------------------------------------------------------
USE green
@@ -161,7 +169,7 @@ PROGRAM relax
IMPLICIT NONE
REAL*8, PARAMETER :: DEG2RAD = 0.01745329251994329547437168059786927_8
- INTEGER, PARAMETER :: ITERATION_MAX = 900
+ INTEGER, PARAMETER :: ITERATION_MAX = 9900
REAL*8, PARAMETER :: STEP_MAX = 1e7
INTEGER :: i,k,sx1,sx2,sx3,e,ne,nv,np,nop,npl,nps,oi,nfc, &
@@ -183,7 +191,7 @@ PROGRAM relax
CHARACTER(80) :: rffilename,vcfilename,cgfilename
CHARACTER(3) :: digit
#endif
- REAL*8 :: dx1,dx2,dx3,oz,ozs,t,Dt,tm,odt
+ REAL*8 :: dx1,dx2,dx3,oz,ozs,t,Dt,tm,odt,tscale
! coseismic events
TYPE(EVENT_STRUC), DIMENSION(:), ALLOCATABLE :: events
TYPE(EVENT_STRUC) :: inter
@@ -371,14 +379,16 @@ PROGRAM relax
CALL tensorfieldadd(sig,tau,sx1,sx2,sx3/2,c1=0._4,c2=-1._4)
CALL stressupdate(u1,u2,u3,lambda,mu,dx1,dx2,dx3,sx1,sx2,sx3/2,sig)
- ! export stress
+ IF (isoutput(skip,t,i,odt,oi,events(e)%time)) THEN
+ ! export stress
#ifdef GRD
- CALL exportstressgrd(sig,sx1,sx2,sx3/2,dx1,dx2,dx3,ozs,x0,y0,wdir,i-1)
+ CALL exportstressgrd(sig,sx1,sx2,sx3/2,dx1,dx2,dx3,ozs,x0,y0,wdir,i-1)
#endif
#ifdef PROJ
- CALL exportstressproj(sig,sx1,sx2,sx3/2,dx1,dx2,dx3,ozs, &
- x0,y0,lon0,lat0,zone,umult,wdir,i-1)
+ CALL exportstressproj(sig,sx1,sx2,sx3/2,dx1,dx2,dx3,ozs, &
+ x0,y0,lon0,lat0,zone,umult,wdir,i-1)
#endif
+ END IF
! initialize large time step
tm=STEP_MAX;
@@ -410,7 +420,7 @@ PROGRAM relax
IF (ALLOCATED(faultcreepstruc)) THEN
DO k=1,np
CALL frictioneigenstress(n(k)%x,n(k)%y,n(k)%z, &
- n(k)%width,n(k)%length,n(k)%strike,n(k)%dip,beta, &
+ n(k)%width,n(k)%length,n(k)%strike,n(k)%dip,n(k)%rake,beta, &
sig,mu,faultcreepstruc,sx1,sx2,sx3/2,dx1,dx2,dx3,moment, &
maxwelltime=maxwell(3))
END DO
@@ -433,7 +443,7 @@ PROGRAM relax
END IF
! choose an integration time step
- CALL integrationstep(tm,Dt,t,oi,odt,events,e,ne)
+ CALL integrationstep(tm,Dt,t,oi,odt,skip,tscale,events,e,ne)
CALL tensorfieldadd(sig,moment,sx1,sx2,sx3/2,c1=0.0_4,c2=1._4)
@@ -480,17 +490,21 @@ PROGRAM relax
! use v1 as placeholders for the afterslip planes
v1=0
DO k=1,np
+ ! one may use optional arguments ...,VEL=v1) to convert
+ ! fault slip to eigenstrain (scalar)
CALL frictioneigenstress(n(k)%x,n(k)%y,n(k)%z, &
- n(k)%width,n(k)%length,n(k)%strike,n(k)%dip,beta, &
- sig,mu,faultcreepstruc,sx1,sx2,sx3/2,dx1,dx2,dx3,moment,VEL=v1)
+ n(k)%width,n(k)%length,n(k)%strike,n(k)%dip,n(k)%rake,beta, &
+ sig,mu,faultcreepstruc,sx1,sx2,sx3/2,dx1,dx2,dx3,moment)
END DO
! update slip history
CALL fieldadd(gamma,v1,sx1+2,sx2,sx3/2,c2=REAL(Dt))
! export strike and dip creep velocity
- CALL exportcreep(np,n,beta,sig,faultcreepstruc, &
- sx1,sx2,sx3/2,dx1,dx2,dx3,x0,y0,wdir,oi)
+ IF (isoutput(skip,t,i,odt,oi,events(e)%time)) THEN
+ CALL exportcreep(np,n,beta,sig,faultcreepstruc, &
+ sx1,sx2,sx3/2,dx1,dx2,dx3,x0,y0,wdir,oi)
+ END IF
END IF
! interseismic loading
@@ -536,6 +550,7 @@ PROGRAM relax
IF (e .LT. ne) THEN
IF (abs(t-events(e+1)%time) .LT. 1e-6) THEN
e=e+1
+ events(e)%i=i
PRINT '("coseismic event ",I3.3)', e
PRINT 0990
@@ -555,7 +570,7 @@ PROGRAM relax
END IF
END IF
- ! points are exported systematically
+ ! points are exported at all time steps
IF (ALLOCATED(ptsname)) THEN
CALL exportpoints(u1,u2,u3,sx1,sx2,sx3/2,dx1,dx2,dx3, &
opts,ptsname,t,wdir,.false.,x0,y0,rot)
@@ -851,7 +866,7 @@ CONTAINS
#endif
INTEGER :: iunit
!$ INTEGER :: omp_get_num_procs,omp_get_max_threads
- REAL*8 :: dummy
+ REAL*8 :: dummy,dum1,dum2
! default is standard input
IF (.NOT. PRESENT(unit)) THEN
@@ -919,7 +934,7 @@ CONTAINS
WRITE_DEBUG_INFO
WRITE (0,'("invalid UTM zone ",I," (1<=zone<=60. exiting.)")') zone
STOP 1
- ENDIF
+ END IF
#endif
PRINT '(a)', "observation depth (displacement and stress)"
@@ -962,14 +977,16 @@ CONTAINS
READ (dataline,*) lambda,mu,gam
PRINT '(3ES10.2E2)',lambda,mu,gam
- PRINT '(a)', "integration time and time step"
+ PRINT '(a)', "time interval, (positive time step) or (negative skip, scaling)"
CALL getdata(unit,dataline)
READ (dataline,*) interval, odt
IF (odt .LT. 0.) THEN
- skip=fix(-odt)
- PRINT '(ES9.2E1," (output every ",I3.3," computational steps)")', interval,skip
+ READ (dataline,*) dum1, dum2, tscale
+ skip=ceiling(-odt)
+ PRINT '(ES9.2E1," (output every ",I3.3," steps, dt scaled by ",ES7.2E1,")")', &
+ interval,skip,tscale
ELSE
- PRINT '(2ES9.2E1)', interval,odt
+ PRINT '(ES9.2E1," (output every ",ES9.2E1," time unit)")', interval,odt
END IF
@@ -1338,17 +1355,17 @@ CONTAINS
IF (iostatus>0) STOP "could not allocate the plane list"
PRINT 2000
- PRINT 2100
+ PRINT 2500
PRINT 2000
DO k=1,np
CALL getdata(unit,dataline)
READ (dataline,*) i,n(k)%x,n(k)%y,n(k)%z,&
- n(k)%length,n(k)%width,n(k)%strike,n(k)%dip
+ n(k)%length,n(k)%width,n(k)%strike,n(k)%dip,n(k)%rake
- PRINT '(I3.3," ",5ES9.2E1,2f7.1)',i, &
+ PRINT '(I3.3," ",5ES9.2E1,3f7.1)',i, &
n(k)%x,n(k)%y,n(k)%z, &
- n(k)%length,n(k)%width,n(k)%strike,n(k)%dip
+ n(k)%length,n(k)%width,n(k)%strike,n(k)%dip,n(k)%rake
IF (i .ne. k) THEN
WRITE_DEBUG_INFO
@@ -1358,7 +1375,7 @@ CONTAINS
! comply to Wang's convention
CALL wangconvention(dummy,n(k)%x,n(k)%y,n(k)%z,&
- n(k)%length,n(k)%width,n(k)%strike,n(k)%dip,dummy,rot)
+ n(k)%length,n(k)%width,n(k)%strike,n(k)%dip,n(k)%rake,rot)
#ifdef VTK
! export the afterslip segment in VTK format
@@ -1535,6 +1552,7 @@ CONTAINS
END IF
ELSE
events(1)%time=0._8
+ events(1)%i=0
END IF
! - - - - - - - - - - - - - - - - - - - - - - - - - -
@@ -1607,6 +1625,8 @@ CONTAINS
rffilename=wdir(1:j-1)//"/rfaults-"//digit//".vtp"
CALL exportvtk_rfaults(events(e),rffilename)
#endif
+ rffilename=wdir(1:j-1)//"/rfaults-"//digit//".dat"
+ CALL exportxy_rfaults(events(e),rffilename)
PRINT 2000
END IF
@@ -1764,6 +1784,7 @@ 2200 FORMAT ("no. slip x1
2200 FORMAT ("no. slip x1 x2 x3 length width strike dip rake")
2300 FORMAT ("no. name x1 x2 x3 (name is a 4-character string)")
2400 FORMAT ("no. strain x1 x2 x3 (positive for extension)")
+2500 FORMAT ("no. x1 x2 x3 length width strike dip rake")
END SUBROUTINE init
@@ -1784,7 +1805,7 @@ 2400 FORMAT ("no. strain x1
INTEGER, INTENT(IN) :: skip,i,oi
REAL*8, INTENT(IN) :: t,odt,etime
- IF (((0 .EQ. skip) .AND. (abs(t-oi*odt) .LT. 1e-6)) .OR. &
+ IF (((0 .EQ. skip) .AND. (abs(t-oi*odt) .LT. 1e-6*odt)) .OR. &
((0 .LT. skip) .AND. (MOD(i-1,skip) .EQ. 0)) .OR. &
(abs(t-etime) .LT. 1e-6)) THEN
isoutput=.TRUE.
@@ -1796,29 +1817,69 @@ 2400 FORMAT ("no. strain x1
!--------------------------------------------------------------------
! subroutine IntegrationStep
- ! find the time-integration forward step based on user-defined
- ! conditions. by default, time step is five times smaller than the
- ! instantaneous Maxwell relaxation time. Time step can be reduced
- ! so that next step corresponds to a following coseismic event.
+ ! find the time-integration forward step for the predictor-corrector
+ ! scheme.
+ !
+ ! input file line
+ !
+ ! time interval, (positive dt step) or (negative skip and scaling)
+ !
+ ! can be filled by either 1)
+ !
+ ! T, dt
+ !
+ ! where T is the time interval of the simulation and dt is the
+ ! output time step, or 2)
+ !
+ ! T, -n, t_s
+ !
+ ! where n indicates the number of computational steps before
+ ! outputing results, t_s is a scaling applied to internally
+ ! computed time step.
+ !
+ ! for case 1), an optimal time step is evaluated internally to
+ ! ensure stability (t_m/10) of time integration. The actual
+ ! time step Dt is chosen as
+ !
+ ! Dt = min( t_m/10, ((t%odt)+1)*odt-t )
+ !
+ ! where t is the current time in the simulation. regardless of
+ ! time step Dt, results are output if t is a multiple of dt.
+ !
+ ! for case 2), the time step is chosen internally based on an
+ ! estimate of the relaxation time (t_m/10). Results are output
+ ! every n steps. The actual time step is chosen as
+ !
+ ! Dt = min( t_m/10*t_s, t(next event)-t )
+ !
+ ! where index is the number of computational steps after a coseismic
+ ! event and t(next event) is the time of the next coseismic event.
!
! sylvain barbot (01/01/08) - original form
!--------------------------------------------------------------------
- SUBROUTINE integrationstep(tm,Dt,t,oi,odt,events,e,ne)
- REAL*8, INTENT(INOUT) :: tm,Dt
- REAL*8, INTENT(IN) :: t,odt
- INTEGER, INTENT(IN) :: oi,e,ne
+ SUBROUTINE integrationstep(tm,Dt,t,oi,odt,skip,tscale,events,e,ne)
+ REAL*8, INTENT(INOUT) :: tm,Dt,odt
+ REAL*8, INTENT(IN) :: t,tscale
+ INTEGER, INTENT(IN) :: oi,e,ne,skip
TYPE(EVENT_STRUC), INTENT(IN), DIMENSION(:) :: events
+ ! output at optimal computational intervals
Dt=tm/10._8
+
+ ! reduce time in case something happens in [ t, t+Dt ]
IF (0 .EQ. skip) THEN
- ! uniform output interval
- IF ((t+Dt) .GE. (dble(oi)*odt)-Dt*0.04) THEN
+ ! reduce time step so that t+Dt is time at next
+ ! user-required output time
+ IF ((t+Dt) .GE. (dble(oi)*odt)-Dt*0.04d0) THEN
! pick a smaller time step to reach :
! integers of odt
Dt=dble(oi)*odt-t
END IF
ELSE
- ! output at optimal computational intervals
+ ! scale the estimate of optimal time step
+ Dt=Dt*tscale
+
+ ! reduce time step so that t+Dt is time to next event
IF (e .LT. ne) THEN
IF ((t+Dt-events(e+1)%time) .GE. 0._8) THEN
! pick a smaller time step to reach
More information about the CIG-COMMITS
mailing list