[cig-commits] commit: export afterslip velocity to .vtk; change the name of afterslip output files.
Mercurial
hg at geodynamics.org
Sat Apr 7 17:03:48 PDT 2012
changeset: 75:914810ee9240
user: Sylvain Barbot <sbarbot at caltech.edu>
date: Sun Apr 01 19:03:08 2012 -0700
files: examples/tutorials/run3.input src/export.f90 src/friction3d.f90 src/relax.f90 util/grdmap.sh
description:
export afterslip velocity to .vtk; change the name of afterslip output files.
diff -r e7295294f654 -r 914810ee9240 examples/tutorials/run3.input
--- a/examples/tutorials/run3.input Sun Apr 01 14:02:51 2012 -0700
+++ b/examples/tutorials/run3.input Sun Apr 01 19:03:08 2012 -0700
@@ -35,16 +35,17 @@ 0
# number of friction interfaces
1
# no depth gamma0 (a-b)sig friction cohesion
- 1 0.3 1e3 1e3 0.0 0
+ 1 0.0 1e3 1e3 0.0 0
# number of creeping faults
-1
+2
# define the fault planes where afterslip occurs. it is possible to
# constrain the rake slip to be within 180 degrees of given direction by
# providing a rake between -180 and 180. for values outside this range,
# the constraint of rake is not applied. this afterslip plane starts
# immediately below the coseismic rupture.
-# no x1 x2 x3 length width strike dip rake (slip rake is +-90 degrees from this value)
- 1 -2 0 0 4 2 0 90 0
+# nb x1 x2 x3 length width strike dip rake (slip rake is +-90 degrees from this value)
+ 1 -2 0.1 0 4 2 0 90 0
+ 2 -2 1.4 0 4 2 0 90 0
# number of interseismic loading strike-slip
0
# number of interseismic loading opening cracks
@@ -53,8 +54,8 @@ 1
1
# number of shear dislocations
1
-# no slip x1 x2 x3 length width strike dip rake
- 1 +1 -0.5 0 0 1 1 0 90 0
+# no slip x1 x2 x3 length width strike dip rake
+ 1 +1 -0.5 0.1 0 1 1 0 90 0
# number of tensile cracks
0
# number of dilatation sources
diff -r e7295294f654 -r 914810ee9240 src/export.f90
--- a/src/export.f90 Sun Apr 01 14:02:51 2012 -0700
+++ b/src/export.f90 Sun Apr 01 19:03:08 2012 -0700
@@ -899,34 +899,90 @@ END SUBROUTINE exporteigenstrain
WRITE (digit,'(I3.3)') i
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,n(k)%rake,beta, &
- sx1,sx2,sx3,dx1,dx2,dx3,sig,structure,n(k)%patch)
ns1=SIZE(n(k)%patch,1)
ns2=SIZE(n(k)%patch,2)
- !patch(:,:)%x1=patch(:,:)%x1+x0
- !patch(:,:)%x2=patch(:,:)%x2+y0
-
WRITE (sdigit,'(I5.5)') k
#ifdef TXT_EXPORTCREEP
- outfile=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".creep.txt"
+ outfile=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".creep.dat"
OPEN (UNIT=15,FILE=outfile,IOSTAT=iostatus,FORM="FORMATTED")
IF (iostatus>0) STOP "could not open file for export"
WRITE (15,'("# x1 x2 x3 yr yz", &
- " slip strike-slip dip-slip")')
- WRITE (15,'(8ES11.3E2)') ((n(k)%patch(i1,i2), i1=1,ns1,skip), i2=1,ns2,skip)
-
+ " slip strike-slip dip-slip velocity ss vel", &
+ " ds vel taus sig11 sig12 sig13 sig22 sig23 sig33")')
+ WRITE (15,'(18ES11.3E2)') ((n(k)%patch(i1,i2)%x1,n(k)%patch(i1,i2)%x3,n(k)%patch(i1,i2)%x3, &
+ n(k)%patch(i1,i2)%lx,n(k)%patch(i1,i2)%lz, &
+ n(k)%patch(i1,i2)%slip, &
+ n(k)%patch(i1,i2)%ss, &
+ n(k)%patch(i1,i2)%ds, &
+ n(k)%patch(i1,i2)%v, &
+ n(k)%patch(i1,i2)%vss, &
+ n(k)%patch(i1,i2)%vds, &
+ n(k)%patch(i1,i2)%taus, &
+ n(k)%patch(i1,i2)%sig%s11, &
+ n(k)%patch(i1,i2)%sig%s12, &
+ n(k)%patch(i1,i2)%sig%s13, &
+ n(k)%patch(i1,i2)%sig%s22, &
+ n(k)%patch(i1,i2)%sig%s23, &
+ n(k)%patch(i1,i2)%sig%s33, &
+ i1=1,ns1,skip), i2=1,ns2,skip)
CLOSE(15)
#endif
+ outfile=wdir(1:pos-1)//"/creep-"//sdigit//"s-"//digit//".vtk"
+
+ OPEN (UNIT=15,FILE=outfile,IOSTAT=iostatus,FORM="FORMATTED")
+ IF (iostatus>0) STOP "could not open file for export"
+
+ WRITE (15,'("# vtk DataFile Version 3.0")')
+ WRITE (15,'("test")')
+ WRITE (15,'("ASCII")')
+ WRITE (15,'("DATASET POLYDATA")')
+ WRITE (15,'("POINTS ",i0," float")') ns1*ns2
+ WRITE (15,'(3ES11.3E2)') ((n(k)%patch(i1,i2)%x1,n(k)%patch(i1,i2)%x2,n(k)%patch(i1,i2)%x3,i1=1,ns1), i2=1,ns2)
+ WRITE (15,'("")')
+ WRITE (15,'("POLYGONS ",i0,X,i0)') (ns1-1)*(ns2-1)*2,(ns1-1)*(ns2-1)*2*4
+ DO i2=1,ns2-1
+ DO i1=1,ns1-1
+ WRITE (15,'("3 ",i0,X,i0,X,i0)') (i2-1)*ns1+i1-1,(i2-1)*ns1+i1-1+1,i2*ns1+i1-1+1
+ WRITE (15,'("3 ",i0,X,i0,X,i0)') (i2-1)*ns1+i1-1,(i2 )*ns1+i1-1 ,i2*ns1+i1-1+1
+ END DO
+ END DO
+ WRITE (15,'("")')
+ WRITE (15,'("POINT_DATA ",i0)') (ns1)*(ns2)
+ WRITE (15,'("SCALARS afterslip float")')
+ WRITE (15,'("LOOKUP_TABLE default")')
+ WRITE (15,'(f)') ((n(k)%patch(i1,i2)%slip,i1=1,ns1),i2=1,ns2)
+ WRITE (15,'("SCALARS strike_slip float")')
+ WRITE (15,'("LOOKUP_TABLE default")')
+ WRITE (15,'(f)') ((n(k)%patch(i1,i2)%ss,i1=1,ns1),i2=1,ns2)
+ WRITE (15,'("SCALARS dip-slip float")')
+ WRITE (15,'("LOOKUP_TABLE default")')
+ WRITE (15,'(f)') ((n(k)%patch(i1,i2)%ds,i1=1,ns1),i2=1,ns2)
+ WRITE (15,'("SCALARS velocity float")')
+ WRITE (15,'("LOOKUP_TABLE default")')
+ WRITE (15,'(f)') ((n(k)%patch(i1,i2)%v,i1=1,ns1),i2=1,ns2)
+ WRITE (15,'("SCALARS velocity_strike_slip float")')
+ WRITE (15,'("LOOKUP_TABLE default")')
+ WRITE (15,'(f)') ((n(k)%patch(i1,i2)%vss,i1=1,ns1),i2=1,ns2)
+ WRITE (15,'("SCALARS velocity_dip_slip float")')
+ WRITE (15,'("LOOKUP_TABLE default")')
+ WRITE (15,'(f)') ((n(k)%patch(i1,i2)%vds,i1=1,ns1),i2=1,ns2)
+ WRITE (15,'("SCALARS shear_stress float")')
+ WRITE (15,'("LOOKUP_TABLE default")')
+ WRITE (15,'(f)') ((n(k)%patch(i1,i2)%taus,i1=1,ns1),i2=1,ns2)
+ CLOSE(15)
+
#ifdef GRD_EXPORTCREEP
- file1=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".creep-north.grd"
- file2=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".creep-east.grd"
- file3=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".creep-up.grd"
+ ALLOCATE(temp1(ns1,ns2),temp2(ns1,ns2),temp3(ns1,ns2),STAT=iostatus)
+ IF (iostatus>0) STOP "could not allocate temporary arrays for GRD slip export."
+
+ file1=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".slip-strike.grd"
+ file2=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".slip-dip.grd"
+ file3=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".slip.grd"
! convert to c standard
j=INDEX(file1," ")
@@ -936,14 +992,44 @@ END SUBROUTINE exporteigenstrain
j=INDEX(file3," ")
file3(j:j)=char(0)
- ALLOCATE(temp1(ns1,ns2),temp2(ns1,ns2),temp3(ns1,ns2),STAT=iostatus)
- IF (iostatus>0) STOP "could not allocate temporary arrays for GRD slip export."
-
DO i2=1,ns2
DO i1=1,ns1
temp1(ns1+1-i1,i2)=n(k)%patch(i1,i2)%ds
temp2(ns1+1-i1,i2)=n(k)%patch(i1,i2)%ss
temp3(ns1+1-i1,i2)=n(k)%patch(i1,i2)%slip
+ END DO
+ END DO
+
+ ! xmin is the lowest coordinates (positive eastward in GMT)
+ xmin= MINVAL(n(k)%patch(:,:)%lx)
+ ! ymin is the lowest coordinates (positive northward in GMT)
+ ymin=-MAXVAL(n(k)%patch(:,:)%lz)
+
+ ! call the c function "writegrd_"
+ CALL writegrd(temp1,ns1,ns2,ymin,xmin,dx3,dx2, &
+ rland,rdum,title,file1)
+ CALL writegrd(temp2,ns1,ns2,ymin,xmin,dx3,dx2, &
+ rland,rdum,title,file2)
+ CALL writegrd(temp3,ns1,ns2,ymin,xmin,dx3,dx2, &
+ rland,rdum,title,file3)
+
+ file1=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".vel-strike.grd"
+ file2=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".vel-dip.grd"
+ file3=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".vel.grd"
+
+ ! convert to c standard
+ j=INDEX(file1," ")
+ file1(j:j)=char(0)
+ j=INDEX(file2," ")
+ file2(j:j)=char(0)
+ j=INDEX(file3," ")
+ file3(j:j)=char(0)
+
+ DO i2=1,ns2
+ DO i1=1,ns1
+ temp1(ns1+1-i1,i2)=n(k)%patch(i1,i2)%vds
+ temp2(ns1+1-i1,i2)=n(k)%patch(i1,i2)%vss
+ temp3(ns1+1-i1,i2)=n(k)%patch(i1,i2)%v
END DO
END DO
@@ -1153,7 +1239,7 @@ END SUBROUTINE exportcreep
WRITE (15,'(" <Piece NumberOfPoints=",a,"8",a," NumberOfPolys=",a,"6",a,">")'),q,q,q,q
WRITE (15,'(" <Points>")')
WRITE (15,'(" <DataArray type=",a,"Float32",a, &
- " Name=",a,"Comp. Grid",a, &
+ " Name=",a,"Computational Grid",a, &
" NumberOfComponents=",a,"3",a, &
" format=",a,"ascii",a,">")'),q,q,q,q,q,q,q,q
WRITE (15,'(24ES9.2E1)') &
diff -r e7295294f654 -r 914810ee9240 src/friction3d.f90
--- a/src/friction3d.f90 Sun Apr 01 14:02:51 2012 -0700
+++ b/src/friction3d.f90 Sun Apr 01 19:03:08 2012 -0700
@@ -454,8 +454,12 @@ CONTAINS
CALL local2ref(xr,yr,zr,x1,x2,x3)
! initialize zero slip velocity
- patch(j2,j3)=SLIPPATCH_STRUCT(x1,x2,x3,yr,zr,0._8,0._8,0._8, &
- 0._8,0._8,0._8,0._8,s)
+ patch(j2,j3)%x1=x1
+ patch(j2,j3)%x2=x2
+ patch(j2,j3)%x3=x3
+ patch(j2,j3)%lx=yr
+ patch(j2,j3)%lz=zr
+ patch(j2,j3)%sig=s
! discard out-of-bound locations
IF ( (x1 .GT. DBLE(sx1/2-1)*dx1) .OR. (x1 .LT. -DBLE(sx1/2)*dx1) &
@@ -492,18 +496,15 @@ CONTAINS
patch(j2,j3)%taus=taus
! creep rate
- patch(j2,j3)%slip=vo*2._8*sinh(tau/tauc)
patch(j2,j3)%v=vo*2._8*sinh(tau/tauc)
! shear traction direction
ts=ts/taus
! strike-direction creep rate
- patch(j2,j3)%ss=patch(j2,j3)%slip*SUM(ts*sv)
patch(j2,j3)%vss=patch(j2,j3)%v*SUM(ts*sv)
! dip-direction creep rate
- patch(j2,j3)%ds=patch(j2,j3)%slip*SUM(ts*dv)
patch(j2,j3)%vds=patch(j2,j3)%v*SUM(ts*dv)
END DO
@@ -550,4 +551,41 @@ CONTAINS
END SUBROUTINE monitorfriction
+ !---------------------------------------------------------------------
+ !> function FrictionAdd
+ !! update the cumulative slip of a creeping segment
+ !!
+ !! \author sylvain barbot (04-01-12) - original form
+ !---------------------------------------------------------------------
+ SUBROUTINE frictionadd(np,n,dt)
+ INTEGER, INTENT(IN) :: np
+ TYPE(PLANE_STRUCT), INTENT(INOUT), DIMENSION(np) :: n
+ REAL*8, INTENT(IN) :: dt
+
+ INTEGER :: px2,px3,j2,j3,k
+
+ ! number of samples in the dip and strike direction
+ DO k=1,np
+ px2=SIZE(n(k)%patch,1)
+ px3=SIZE(n(k)%patch,2)
+
+ ! loop in the dip direction
+ DO j3=1,px3
+ ! loop in the strike direction
+ DO j2=1,px2
+ ! cumulative creep
+ n(k)%patch(j2,j3)%slip=n(k)%patch(j2,j3)%slip+dt*n(k)%patch(j2,j3)%v
+
+ ! cumulative strike-direction creep
+ n(k)%patch(j2,j3)%ss=n(k)%patch(j2,j3)%ss+dt*n(k)%patch(j2,j3)%vss
+
+ ! cumulative dip-direction creep
+ n(k)%patch(j2,j3)%ds=n(k)%patch(j2,j3)%ds+dt*n(k)%patch(j2,j3)%vds
+ END DO
+ END DO
+
+ END DO
+
+ END SUBROUTINE frictionadd
+
END MODULE friction3d
diff -r e7295294f654 -r 914810ee9240 src/relax.f90
--- a/src/relax.f90 Sun Apr 01 14:02:51 2012 -0700
+++ b/src/relax.f90 Sun Apr 01 19:03:08 2012 -0700
@@ -606,13 +606,14 @@ PROGRAM relax
in%n(k)%strike,in%n(k)%dip,in%n(k)%rake,in%beta, &
sig,in%mu,in%faultcreepstruc,in%sx1,in%sx2,in%sx3/2, &
in%dx1,in%dx2,in%dx3,moment)
+
+ ! keep track if afterslip instantaneous velocity
+ CALL monitorfriction(in%n(k)%x,in%n(k)%y,in%n(k)%z, &
+ in%n(k)%width,in%n(k)%length, &
+ in%n(k)%strike,in%n(k)%dip,in%n(k)%rake,in%beta, &
+ in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
+ sig,in%faultcreepstruc,in%n(k)%patch)
END DO
-
- ! export strike and dip creep velocity
- IF (isoutput(in%skip,t,i,in%odt,oi,in%events(e)%time)) THEN
- CALL exportcreep(in%np,in%n,in%beta,sig,in%faultcreepstruc, &
- in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%x0,in%y0,in%wdir,oi)
- END IF
END IF
@@ -658,6 +659,7 @@ PROGRAM relax
CALL fieldadd(u2,v2,in%sx1+2,in%sx2,in%sx3/2,c2=REAL(Dt))
CALL fieldadd(u3,v3,in%sx1+2,in%sx2,in%sx3/2,c2=REAL(Dt))
CALL tensorfieldadd(tau,moment,in%sx1,in%sx2,in%sx3/2,c2=REAL(Dt))
+ CALL frictionadd(in%np,in%n,Dt)
! keep track of the viscoelastic contribution alone
IF (in%isoutputrelax) THEN
@@ -721,6 +723,10 @@ PROGRAM relax
CALL reporttime(1,t,in%reporttimefilename)
+ ! export strike and dip creep velocity
+ CALL exportcreep(in%np,in%n,in%beta,sig,in%faultcreepstruc, &
+ in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%x0,in%y0,in%wdir,oi)
+
! export
#ifdef TXT
IF (in%isoutputtxt) THEN
@@ -735,7 +741,7 @@ PROGRAM relax
END IF
END IF
#endif
- CALL exporteigenstrain(gamma,in%nop,in%op,in%x0,in%y0,in%dx1,in%dx2,in%dx3,in%sx1,in%sx2,in%sx3/2,in%wdir,oi)
+ !CALL exporteigenstrain(gamma,in%nop,in%op,in%x0,in%y0,in%dx1,in%dx2,in%dx3,in%sx1,in%sx2,in%sx3/2,in%wdir,oi)
#ifdef GRD
IF (in%isoutputgrd) THEN
IF (in%isoutputrelax) THEN
diff -r e7295294f654 -r 914810ee9240 util/grdmap.sh
--- a/util/grdmap.sh Sun Apr 01 14:02:51 2012 -0700
+++ b/util/grdmap.sh Sun Apr 01 19:03:08 2012 -0700
@@ -68,16 +68,55 @@ rm -f $colors
}
-if (test $# -lt "1"); then
- echo "usage: grdmap.sh [-b -5/5/-5/5] [-p -1.5/1.5/0.1] [...] file1.grd ... fileN.grd"
+gmtset LABEL_FONT_SIZE 12p
+gmtset HEADER_FONT_SIZE 12p
+gmtset ANOT_FONT_SIZE 12p
+gmtset COLOR_BACKGROUND 0/0/255
+gmtset COLOR_FOREGROUND 255/0/0
+gmtset COLOR_NAN 255/255/255
+gmtset PAPER_MEDIA a5
+
+libdir=$(dirname $0)
+EXTRA=""
+
+while getopts "b:c:e:ghi:o:p:v:s:t:u:xrY:" flag
+do
+ case "$flag" in
+ b) bset=1;bds=$OPTARG;;
+ c) cset=1;carg=$OPTARG;;
+ e) eset=1;EXTRA="$EXTRA $OPTARG";;
+ g) gset=1;;
+ h) hset=1;;
+ i) iset=1;illumination="-I$OPTARG";;
+ o) oset=1;ODIR=$OPTARG;;
+ p) pset=1;P=$OPTARG;PSSCALE=`echo $OPTARG | awk -F"/" 'function abs(x){return x<0?-x:x}{print abs($2-$1)/4}'`;echo $PSSCALE;;
+ v) vset=1;SIZE=$OPTARG;VECTOR=$OPTARG"c";;
+ r) rset=1;;
+ s) sset=1;ADX=$OPTARG;;
+ t) tset=1;tick=$OPTARG;;
+ u) uset=1;unit=$OPTARG;;
+ x) xset=1;;
+ Y) Yset=1;Yshift=$OPTARG;;
+ esac
+done
+for item in $bset $cset $iset $oset $pset $vset $sset $tset $uset $Yset $EXTRA;do
+ shift;shift
+done
+for item in $gset $hset $xset $rset;do
+ shift;
+done
+
+if [ $# -lt "1" ] || [ "$hset" == "1" ] ; then
+ echo "usage: map.sh [-b -5/5/-5/5] [-p -1.5/1.5/0.1] [...] file1.grd ... fileN.grd"
echo "or"
- echo " grdmap.sh [-b -5/5/-5/5] [-p -1.5/1.5/0.1] [...] dir1/index1 ... dirN/indexN"
+ echo " map.sh [-b -5/5/-5/5] [-p -1.5/1.5/0.1] [...] dir1/index1 ... dirN/indexN"
echo ""
echo "options:"
echo " -b bds sets the map bound to bds"
echo " -c palette_name [default my_jet]"
echo " -e file.sh runs file.sh to add content to map"
echo " -g switch to geographic projection (longitude/latitude)"
+ echo " -h display this error message and exit"
echo " -i file.grd illuminates the grd image with file.grd"
echo " -o output directory"
echo " -p palette sets the color scale bounds to palette"
@@ -95,43 +134,6 @@ if (test $# -lt "1"); then
exit
fi
-
-gmtset LABEL_FONT_SIZE 12p
-gmtset HEADER_FONT_SIZE 12p
-gmtset ANOT_FONT_SIZE 12p
-gmtset COLOR_BACKGROUND 0/0/255
-gmtset COLOR_FOREGROUND 255/0/0
-gmtset COLOR_NAN 255/255/255
-gmtset PAPER_MEDIA a5
-
-libdir=$(dirname $0)
-EXTRA=""
-
-while getopts "b:c:e:gi:o:p:v:s:t:u:xrY:" flag
-do
- case "$flag" in
- b) bset=1;bds=$OPTARG;;
- c) cset=1;carg=$OPTARG;;
- e) eset=1;EXTRA="$EXTRA $OPTARG";;
- g) gset=1;;
- i) iset=1;illumination="-I$OPTARG";;
- o) oset=1;ODIR=$OPTARG;;
- p) pset=1;P=$OPTARG;PSSCALE=`echo $OPTARG | awk -F"/" 'function abs(x){return x<0?-x:x}{print abs($2-$1)/4}'`;echo $PSSCALE;;
- v) vset=1;SIZE=$OPTARG;VECTOR=$OPTARG"c";;
- r) rset=1;;
- s) sset=1;ADX=$OPTARG;;
- t) tset=1;tick=$OPTARG;;
- u) uset=1;unit=$OPTARG;;
- x) xset=1;;
- Y) Yset=1;Yshift=$OPTARG;;
- esac
-done
-for item in $bset $cset $iset $oset $pset $vset $sset $tset $uset $Yset $EXTRA;do
- shift;shift
-done
-for item in $gset $xset $rset;do
- shift;
-done
while [ "$#" != "0" ];do
More information about the CIG-COMMITS
mailing list