[cig-commits] commit: output time series of Coulomb stress.
Mercurial
hg at geodynamics.org
Wed Oct 19 15:18:58 PDT 2011
changeset: 31:8077fdb82342
tag: tip
user: Sylvain Barbot <sbarbot at caltech.edu>
date: Wed Oct 05 20:52:55 2011 -0700
files: examples/run4.sh export.f90 input.f90 readme.txt relax.f90
description:
output time series of Coulomb stress.
diff -r fc4919855571 -r 8077fdb82342 examples/run4.sh
--- a/examples/run4.sh Wed Oct 05 19:18:00 2011 -0700
+++ b/examples/run4.sh Wed Oct 05 20:52:55 2011 -0700
@@ -43,15 +43,13 @@ if [ ! -e $WDIR ]; then
mkdir $WDIR
fi
-time ../relax --no-stress-output $* <<EOF
+time ../relax --no-stress-output --no-proj-output $* <<EOF
# grid size (sx1,sx2,sx3)
256 256 256
# sampling size, smoothing & nyquist (dx1,dx2,dx3,beta,nq)
0.05 0.05 0.05 0.2 1
# origin position & rotation
0 0 0
-# geographic origin (longitude, latitude, UTM zone, unit)
--120 34 11 1e3
# observation depth (displacement and stress) (stress in only exported in GRD)
0 0
# output directory
@@ -65,7 +63,13 @@ 0
# number of observation points
0
# number of Coulomb patches
-0
+5
+# no. x1 x2 x3 length width strike dip friction
+ 1 -1.0 1 0.5 0.1 0.1 0 90 0.6
+ 2 -0.5 1 0.5 0.1 0.1 90 90 0.6
+ 3 0.0 1 0.5 0.1 0.1 0 90 0.6
+ 4 +0.5 1 0.5 0.1 0.1 45 90 0.6
+ 5 +1.0 1 0.5 0.1 0.1 0 90 0.6
# number of prestress interfaces
0
# number of linear viscous interfaces
diff -r fc4919855571 -r 8077fdb82342 export.f90
--- a/export.f90 Wed Oct 05 19:18:00 2011 -0700
+++ b/export.f90 Wed Oct 05 20:52:55 2011 -0700
@@ -1486,6 +1486,13 @@ END SUBROUTINE exportcreep
WRITE (15,'(ES11.2)'), taun
WRITE (15,'(" </DataArray>")')
+ WRITE (15,'(" <DataArray type=",a,"Float32",a, &
+ " Name=",a,"Coulomb stress",a, &
+ " NumberOfComponents=",a,"1",a, &
+ " format=",a,"ascii",a,">")'), q,q,q,q,q,q,q,q
+ WRITE (15,'(ES11.2)'), taucoulomb
+ WRITE (15,'(" </DataArray>")')
+
WRITE (15,'(" </CellData>")')
WRITE (15,'(" </Piece>")')
@@ -1498,6 +1505,140 @@ END SUBROUTINE exportcreep
CLOSE(15)
END SUBROUTINE exportvtk_rfaults_stress
+
+ !--------------------------------------------------------------------------------
+ !> subroutine ExportCoulombStress
+ !! sample the stress tensor, shear and normal stress and Coulomb
+ !! stress at a series of locations.
+ !!
+ !! each fault patch is attributed to a file in which the time
+ !! evolution is listed in the following format:
+ !!
+ !! #t s11 s12 s13 s22 s23 s33 taus taud tau taun Coulomb
+ !! t0 s11(t0) s12(t0) s13(t0) s22(t0) s23(t0) s33(t0) taus(t0) taud(t0) tau(t0) taun(t0) Coulomb(t0)
+ !! t1 s11(t1) s12(t1) s13(t1) s22(t1) s23(t1) s33(t1) taus(t1) taud(t1) tau(t1) taun(t1) Coulomb(t0)
+ !! ...
+ !!
+ !! where sij(t0) is the component ij of the stress tensor at time t0, taus is
+ !! the component of shear in the strike direction, taud is the component of shear
+ !! in the fault dip direction, tau^2=taus^2+taud^2, taun is the fault normal
+ !! stress and Coulomb(t0) is the Coulomb stress tau+mu*taun.
+ !!
+ !! \author sylvain barbot (10/05/11) - original form
+ !--------------------------------------------------------------------------------
+ SUBROUTINE exportcoulombstress(sig,sx1,sx2,sx3,dx1,dx2,dx3, &
+ nsop,sop,time,wdir,isnew)
+ INTEGER, INTENT(IN) :: sx1,sx2,sx3,nsop
+ TYPE(TENSOR), INTENT(IN), DIMENSION(sx1,sx2,sx3) :: sig
+ TYPE(SEGMENT_STRUCT), INTENT(INOUT), DIMENSION(nsop) :: sop
+ REAL*8, INTENT(IN) :: dx1,dx2,dx3,time
+ CHARACTER(80), INTENT(IN) :: wdir
+ LOGICAL, INTENT(IN) :: isnew
+
+ INTEGER :: iostatus,k,i1,i2,i3
+ CHARACTER :: q
+ CHARACTER(4) :: digit4
+ CHARACTER(80) :: file
+ REAL*8 :: strike,dip,x1,x2,x3,cstrike,sstrike,cdip,sdip,L,W
+ ! segment normal vector, strike direction, dip direction
+ REAL*8, DIMENSION(3) :: n,s,d
+ ! local value of stress
+ TYPE(TENSOR) :: lsig
+ ! stress components
+ REAL*8 :: taun,taus,taustrike,taudip,taucoulomb
+ ! friction coefficient
+ REAL*8 :: friction
+ ! traction components
+ REAL*8, DIMENSION(3) :: t,ts
+
+ IF (0.GE.nsop) RETURN
+
+ ! double-quote character
+ q=char(34)
+
+ DO k=1,nsop
+ WRITE (digit4,'(I4.4)') k
+ file=trim(wdir)//"/cfaults-sigma-"//digit4//".txt"
+
+ ! fault center position
+ x1=sop(k)%x
+ x2=sop(k)%y
+ x3=sop(k)%z
+
+ IF (isnew) THEN
+ OPEN (UNIT=15,FILE=file,IOSTAT=iostatus,FORM="FORMATTED")
+ WRITE (15,'("# center position (north, east, down): ",3ES9.2)') x1,x2,x3
+ WRITE (15,'("# t s11 s12 s13 ", &
+ "s22 s23 s33 taus taud tau taun Coulomb")')
+ ELSE
+ OPEN (UNIT=15,FILE=file,POSITION="APPEND",&
+ IOSTAT=iostatus,FORM="FORMATTED")
+ END IF
+ IF (iostatus>0) STOP "could not open point file for writing"
+
+ ! friction coefficient
+ friction=sop(k)%friction
+
+ ! fault orientation
+ strike=sop(k)%strike
+ dip=sop(k)%dip
+
+ CALL shiftedindex(x1,x2,x3,sx1,sx2,sx3,dx1,dx2,dx3,i1,i2,i3)
+ lsig=sig(i1,i2,i3)
+
+ ! fault dimension
+ W=sop(k)%width
+ L=sop(k)%length
+
+ cstrike=cos(strike)
+ sstrike=sin(strike)
+ cdip=cos(dip)
+ sdip=sin(dip)
+
+ ! surface normal vector components
+ n(1)=+cdip*cstrike
+ n(2)=-cdip*sstrike
+ n(3)=-sdip
+
+ ! 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
+
+ ! traction vector
+ t=lsig .tdot. n
+
+ ! signed normal component
+ taun=SUM(t*n)
+
+ ! shear traction
+ ts=t-taun*n
+
+ ! absolute value of shear component
+ taus=SQRT(SUM(ts*ts))
+
+ ! strike-direction shear component
+ taustrike=SUM(ts*s)
+
+ ! dip-direction shear component
+ taudip=SUM(ts*d)
+
+ ! Coulomb stress
+ taucoulomb=taus+friction*taun
+
+ WRITE (15,'(12ES11.3E2)') time, &
+ lsig%s11,lsig%s12,lsig%s13, &
+ lsig%s22,lsig%s23,lsig%s33, &
+ taustrike,taudip,taus,taun,taucoulomb
+ CLOSE(15)
+ END DO
+
+ END SUBROUTINE exportcoulombstress
!------------------------------------------------------------------
!> subroutine ExportVTK_Rectangle
diff -r fc4919855571 -r 8077fdb82342 input.f90
--- a/input.f90 Wed Oct 05 19:18:00 2011 -0700
+++ b/input.f90 Wed Oct 05 20:52:55 2011 -0700
@@ -286,7 +286,7 @@ CONTAINS
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- ! O B S E R V A T I O N P L A N E S ( C O U L O M B S T R E S S)
+ ! O B S E R V A T I O N P L A N E S
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
PRINT '(a)', "number of observation planes"
CALL getdata(iunit,dataline)
diff -r fc4919855571 -r 8077fdb82342 readme.txt
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/readme.txt Wed Oct 05 20:52:55 2011 -0700
@@ -0,0 +1,94 @@
+
+RELAX - time-dependent postseismic deformation with afterslip and viscoelastic flow.
+
+INSTALLATION:
+
+The code is written in Fortran90 and is optimized for the INTEL ifort compiler.
+The gmt 4.5+ library is required to export results to GMT for post-processing.
+
+More information for compilation is available in the makefile.
+
+RUN:
+
+Some examples are available in the examples directory. Look up the *.sh files for
+comments and explanations.
+
+VISUALIZATION:
+
+Many outputs are exported in the General Mapping Tools (GMT) format, deformation
+maps can be obtained with typical GMT post-processing. Download post-processing
+and visualization tools util.tar available online.
+Make sure to have the .ps file viewer gv or the .pdf file viewer xpdf installed.
+Simulations can be visualized in 3D with the free software Paraview (paraview.org).
+
+PURPOSE:
+
+The program RELAX computes nonlinear time-dependent viscoelastic deformation with
+powerlaw rheology and rate-strengthening friction in a cubic, periodic grid due
+to coseismic stress changes, initial stress, surface loads, and/or moving faults.
+
+DESCRIPTION:
+
+Computation is done semi-analytically inside a cartesian grid. The grid is defined
+by its size sx1*sx2*sx3 and the sampling intervals dx1, dx2 and dx3. Rule of thumb
+is to allow for at least five samples per fault length or width, and to have the
+tip of any fault at least 10 fault widths away from any edge of the computational
+grid.
+
+Coseismic stress changes and initial coseismic deformation results from the
+presence of dislocations in the brittle layer. Fault geometry is prescribed
+following Okada or Wang's convention, with the usual slip, strike, dip and rake and
+is converted to a double-couple equivalent body-force analytically. Current
+implementation allows shear fault (strike slip and dip slip), dykes, Mogi source,
+and surface traction. Faults and dykes can be of arbitrary orientation in the half
+space.
+
+
+INPUT:
+
+Static dislocation sources are discretized into a series of planar segments. Slip
+patches are defined in terms of position, orientation, and slip, as illustrated in
+the following figure:
+
+ N (x1)
+ /
+ /| Strike
+ x1,x2,x3 ->@------------------------ (x2)
+ |\ p . \ W
+ :-\ i . \ i
+ | \ l . \ d
+ :90 \ S . \ t
+ |-Dip\ . \ h
+ : \. | Rake \
+ | -------------------------
+ : L e n g t h
+ Z (x3)
+
+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.
+
+REFERENCE:
+
+More information about parameters and constitutive laws can be found in
+
+S. Barbot and Fialko Y., "Fourier-Domain Green Function for an Elastic Semi-
+Infinite Solid under Gravity, with Applications to Earthquake and Volcano
+Deformation", Geophysical Journal International, v. 182, no. 2, pp. 568-582, 2010,
+doi:10.1111/j.1365-246X.2010.04655.x
+
+and
+
+S. Barbot and Fialko Y., "A Unified Continuum Representation of Postseismic
+Relaxation Mechanisms: Semi-Analytic Models of Afterslip, Poroelastic Rebound and
+Viscoelastic Flow", Geophysical Journal International, v. 182, 3, p. 1124-1140,
+2010, doi:10.1111/j.1365-246X.2010.04678.x
+
+Please cite these papers in publications or public presentations when referring
+to this method.
+
+
diff -r fc4919855571 -r 8077fdb82342 relax.f90
--- a/relax.f90 Wed Oct 05 19:18:00 2011 -0700
+++ b/relax.f90 Wed Oct 05 20:52:55 2011 -0700
@@ -406,6 +406,8 @@ PROGRAM relax
CALL exportvtk_rfaults_stress(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
in%nsop,in%sop,filename)
#endif
+ CALL exportcoulombstress(sig,in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
+ in%nsop,in%sop,0._8,in%wdir,.TRUE.)
CALL reporttime(0,0._8,in%reporttimefilename)
PRINT 1101,0,0._8,0._8,0._8,0._8,0._8,in%interval,0._8,tensoramplitude(tau,in%dx1,in%dx2,in%dx3)
@@ -654,7 +656,7 @@ PROGRAM relax
! points are exported at all time steps
IF (ALLOCATED(in%ptsname)) THEN
CALL exportpoints(u1,u2,u3,sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
- in%opts,in%ptsname,t,in%wdir,.false.,in%x0,in%y0,in%rot)
+ in%opts,in%ptsname,t,in%wdir,.FALSE.,in%x0,in%y0,in%rot)
END IF
! output only at discrete intervals (skip=0, odt>0),
@@ -753,6 +755,8 @@ PROGRAM relax
CALL exportvtk_rfaults_stress(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
in%nsop,in%sop,filename,convention=1,sig=sig)
#endif
+ CALL exportcoulombstress(sig,in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
+ in%nsop,in%sop,t,in%wdir,.FALSE.)
PRINT 1101,i,Dt,maxwell,t,in%interval, &
tensoramplitude(moment,in%dx1,in%dx2,in%dx3), &
More information about the CIG-COMMITS
mailing list