[cig-commits] commit: improves grdmap.sh (better color scale, plots zero fields), add option to output the relaxation component in VTK format, fixes the stress output in VTK format, fixes afterslip example (run3.sh) and add comments about optimal time steps, simplify the monitorfriction subroutine, fixes makefile for fast compilation and debugging.
Mercurial
hg at geodynamics.org
Tue Dec 27 16:04:35 PST 2011
changeset: 53:ed3436b7ac71
tag: tip
user: Sylvain Barbot <sbarbot at caltech.edu>
date: Tue Dec 27 16:04:21 2011 -0800
files: examples/run3.sh export.f90 friction3d.f90 input.f90 makefile_imkl relax.f90 types.f90 util/grdmap.sh writevtk.c
description:
improves grdmap.sh (better color scale, plots zero fields), add option to output the relaxation component in VTK format, fixes the stress output in VTK format, fixes afterslip example (run3.sh) and add comments about optimal time steps, simplify the monitorfriction subroutine, fixes makefile for fast compilation and debugging.
diff -r 8d085e25c7e7 -r ed3436b7ac71 examples/run3.sh
--- a/examples/run3.sh Wed Dec 14 15:04:59 2011 -0800
+++ b/examples/run3.sh Tue Dec 27 16:04:21 2011 -0800
@@ -4,7 +4,8 @@
#
# afterslip following left-lateral strike-slip on a vertical fault.
#
-# output at each computation step (Dt = -1).
+# output at each computation step (Dt = -1) with a reduced time step (1/100 factor) for
+# better numerical accuracy.
#
# to visualize coseismic deformation (requires GRD output, or manual conversion to GRD format):
#
@@ -17,15 +18,66 @@
#
# the -s 0.75 option optimizes the spacing between horizontal vectors.
# to convert .xyz output in .grd output in post processing, type (requires GMT installed on your machine)
-# xyz2grd.sh output3/000
+#
+# xyz2grd.sh output3/000
#
# to visualize a time series of postseismic deformation (requires GRD output):
-# map.sh -b -5/5/-5/5 -p -0.002/0.002/0.0001 -v 0.005 output3/0{01,02,03,04,05,06,07,08,09,10}-relax
+#
+# map.sh -b -5/5/-5/5 -p -0.001/0.001/0.00001 -v 0.001 -s 1 -e rpatch.sh output3/0{00,01,02,03,04,05,06,07,08,09,10}-relax
+#
+# to visualize a time series of fault creep:
+#
+# map.sh -t 1 -s 0.1 -p -0.2/0.2/0.01 output3/0??.s00001.creep-east.grd
#
# type map.sh for a description of command-line options.
#
-# the conversion to geographic is aborted (--no-proj-output) and
+# the conversion to geographic is cancelled (--no-proj-output) and
# the stress component are not exported (--no-stress-output)
+#
+# TIP:
+# the following line in the input file
+#
+# # integration time (t1)
+# 0.5 -1 0.01
+#
+# reduces the time step suggested by the code by a factor of 1/100. This choice
+# of time step leads to slowly decreasing power density
+# |
+# \./
+# I | Dt | tm(ve) | tm(pl) | tm(as) | t/tmax | power | C:E^i |
+# 000* 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00/5.00E-1 0.00E+00 1.03E+00
+# 001* 2.00E-02 1.00E+07 1.00E+07 2.00E+01 2.00E-02/5.00E-1 6.47E-01 1.04E+00
+# 002* 2.00E-02 1.00E+07 1.00E+07 2.00E+01 4.00E-02/5.00E-1 6.11E-01 1.05E+00
+# 003* 2.00E-02 1.00E+07 1.00E+07 2.00E+01 6.00E-02/5.00E-1 5.77E-01 1.06E+00
+# 004* 2.00E-02 1.00E+07 1.00E+07 2.00E+01 8.00E-02/5.00E-1 5.47E-01 1.07E+00
+# 005* 2.00E-02 1.00E+07 1.00E+07 2.00E+01 1.00E-01/5.00E-1 5.18E-01 1.08E+00
+# /^\
+# |
+#
+# It is informative to run the same code with a larger time step, for example
+# with the following parameters:
+#
+# # integration time (t1)
+# 5 -1 0.1
+#
+# which does not lead to numerical convergence.
+# |
+# \./
+# I | Dt | tm(ve) | tm(pl) | tm(as) | t/tmax | power | C:E^i |
+# 000* 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00/5.00E+0 0.00E+00 1.03E+00
+# 001* 2.00E-01 1.00E+07 1.00E+07 2.00E+01 2.00E-01/5.00E+0 4.69E-01 1.12E+00
+# 002* 2.00E-01 1.00E+07 1.00E+07 2.00E+01 4.00E-01/5.00E+0 3.13E-01 1.18E+00
+# 003* 2.00E-01 1.00E+07 1.00E+07 2.00E+01 6.00E-01/5.00E+0 2.18E-01 1.22E+00
+# 004* 2.00E-01 1.00E+07 1.00E+07 2.00E+01 8.00E-01/5.00E+0 1.55E-01 1.26E+00
+# 005* 2.00E-01 1.00E+07 1.00E+07 2.00E+01 1.00E+00/5.00E+0 1.18E-01 1.28E+00
+# 006* 2.00E-01 1.00E+07 1.00E+07 2.00E+01 1.20E+00/5.00E+0 9.96E-02 1.30E+00
+# 007* 2.00E-01 1.00E+07 1.00E+07 2.00E+01 1.40E+00/5.00E+0 9.47E-02 1.31E+00
+# 008* 2.00E-01 1.00E+07 1.00E+07 2.00E+01 1.60E+00/5.00E+0 1.02E-01 1.33E+00
+# 009* 2.00E-01 1.00E+07 1.00E+07 2.00E+01 1.80E+00/5.00E+0 1.35E-01 1.35E+00
+# /^\
+# |
+# in general, a slowly decreasing power density (about 10-20% per step) is
+# a sign of an accurate computation.
WDIR=./output3
@@ -54,9 +106,11 @@ 0 0
# elastic parameters and gamma = (1-nu) rho g / mu = 8.33e-7 /m = 8.33e-4 /km
1 1 8.33e-4
# integration time (t1)
-2 -1 1
+0.5 -1 0.01
# number of observation planes
-0
+1
+# no x1 x2 x3 length width strike dip
+ 1 -3 0 0.5 6 2.5 0 90
# number of observation points
0
# number of Coulomb patches
@@ -70,7 +124,7 @@ 0
# number of friction interfaces
1
# no depth gamma0 (a-b)sig friction cohesion
- 1 1 1e3 1e3 0.6 0
+ 1 1 1e3 1e3 0.0 0
# number of creeping faults
1
# define the fault planes where afterslip occurs. it is possible to
@@ -78,8 +132,8 @@ 1
# 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 +-180 degrees from this value)
- 1 -1 0 1 2 1 0 90 180
+# no x1 x2 x3 length width strike dip rake (slip rake is +-90 degrees from this value)
+ 1 -2 0 1 4 2 0 90 0
# number of interseismic loading strike-slip
0
# number of interseismic loading opening cracks
@@ -88,8 +142,8 @@ 1
1
# number of shear dislocations
1
-# no slip x1 x2 x3 length width strike dip rake
- 1 +1 -1 0 0 2 1 0 90 0
+# no slip x1 x2 x3 length width strike dip rake
+ 1 +1 -0.5 0 0 1 1 0 90 0
# number of tensile cracks
0
# number of dilatation sources
diff -r 8d085e25c7e7 -r ed3436b7ac71 export.f90
--- a/export.f90 Wed Dec 14 15:04:59 2011 -0800
+++ b/export.f90 Tue Dec 27 16:04:21 2011 -0800
@@ -625,14 +625,14 @@ CONTAINS
INTEGER :: iostatus,i1,i2
CHARACTER(80) :: outfiletxt
#endif
-#ifdef GRD_EXPORTEIGENSTRAIN
+!#_indef GRD_EXPORTEIGENSTRAIN
CHARACTER(80) :: outfilegrd
INTEGER :: j,iostat,j1,j2
REAL*4, ALLOCATABLE, DIMENSION(:,:) :: temp
REAL*8 :: rland=9998.,rdum=9999.
REAL*8 :: xmin,ymin
CHARACTER(80) :: title="monitor field "
-#endif
+!#_endif
IF (nop .le. 0) RETURN
@@ -668,7 +668,7 @@ CONTAINS
CLOSE(15)
#endif
-#ifdef GRD_EXPORTEIGENSTRAIN
+!#_ifdef GRD_EXPORTEIGENSTRAIN
outfilegrd=wdir(1:pos-1)//"/"//digit//".s"//sdigit//".estrain.grd"
! convert to c standard
@@ -695,7 +695,7 @@ CONTAINS
DEALLOCATE(temp)
-#endif
+!#_endif
DEALLOCATE(slippatch)
END DO
@@ -741,6 +741,7 @@ END SUBROUTINE exporteigenstrain
!! \author sylvain barbot (01/01/07) - original form
!! (02/25/10) - output in TXT and GRD formats
!---------------------------------------------------------------------
+#define TXT_EXPORTCREEP
SUBROUTINE exportcreep(np,n,beta,sig,structure, &
sx1,sx2,sx3,dx1,dx2,dx3,x0,y0,wdir,i)
INTEGER, INTENT(IN) :: np,sx1,sx2,sx3,i
@@ -790,6 +791,8 @@ END SUBROUTINE exporteigenstrain
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)') ((slippatch(i1,i2), i1=1,ns1,skip), i2=1,ns2,skip)
CLOSE(15)
@@ -1891,7 +1894,6 @@ END SUBROUTINE exportcreep
CHARACTER(80), INTENT(IN) :: filename
INTEGER :: iostatus
- CHARACTER :: q
REAL*8 :: cstrike,sstrike,cdip,sdip
REAL*8, DIMENSION(3) :: s,d,n
diff -r 8d085e25c7e7 -r ed3436b7ac71 friction3d.f90
--- a/friction3d.f90 Wed Dec 14 15:04:59 2011 -0800
+++ b/friction3d.f90 Tue Dec 27 16:04:21 2011 -0800
@@ -283,7 +283,7 @@ CONTAINS
! rake vector component
r(1)=sstrike*cr+cstrike*sdip*sr
r(2)=cstrike*cr-sstrike*sdip*sr
- r(3)=+cdip*sr
+ r(3)=cdip*sr
DO i3=1,sx3
x3=DBLE(i3-1)*dx3
@@ -331,9 +331,9 @@ CONTAINS
! absolute value of shear component
ts=t-taun*n
taus=SQRT(SUM(ts*ts))
-
+
! effective shear stress on fault plane
- tau=taus+friction*taun-cohesion
+ tau=MAX(0.d0,taus+friction*taun-cohesion)
! rake direction test only if | rake | < 3*Pi
IF (SUM(ts*r).LT.0.d0 .AND. ABS(rake).LT.pi2*1.5d0) CYCLE
@@ -351,9 +351,6 @@ CONTAINS
WRITE (0,'("------------------------------------------")')
STOP 5
END IF
-
- ! yield surface test
- IF ((0._8 .GE. taus) .OR. (tau .LE. 0._8)) CYCLE
! shear traction direction
ts=ts/taus
@@ -409,7 +406,7 @@ CONTAINS
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,cr,sr
+ REAL*8 :: cstrike,sstrike,cdip,sdip,cr,sr
REAL*8 :: vo,tauc,taun,taus, &
friction,tau,cohesion
REAL*8 :: x1,x2,x3,xr,yr,zr,Wp,Lp
@@ -449,7 +446,7 @@ CONTAINS
! rake vector component
r(1)=sstrike*cr+cstrike*sdip*sr
r(2)=cstrike*cr-sstrike*sdip*sr
- r(3)=+cdip*sr
+ r(3)=cdip*sr
! loop in the dip direction
DO j3=1,px3+1
@@ -464,60 +461,51 @@ 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)
+
! discard out-of-bound locations
IF ( (x1 .GT. DBLE(sx1/2-1)*dx1) .OR. (x1 .LT. -DBLE(sx1/2)*dx1) &
.OR. (x2 .GT. DBLE(sx2/2-1)*dx2) .OR. (x2 .LT. -DBLE(sx2/2)*dx2) &
- .OR. (x3 .GT. DBLE(sx3-1)*dx3) .OR. (x3 .LT. 0._8) ) THEN
- slip=0._8
- ss=0._8
- ds=0._8
- ELSE
- ! evaluates instantaneous creep velocity
- CALL shiftedindex(x1,x2,x3,sx1,sx2,sx3,dx1,dx2,dx3,i1,i2,i3)
+ .OR. (x3 .GT. DBLE(sx3-1)*dx3) .OR. (x3 .LT. 0._8) ) CYCLE
- ! retrieve friction parameters
- vo=structure(i3)%gammadot0
- tauc=structure(i3)%stressexponent
- friction=structure(i3)%friction
- cohesion=structure(i3)%cohesion
+ ! evaluates instantaneous creep velocity
+ CALL shiftedindex(x1,x2,x3,sx1,sx2,sx3,dx1,dx2,dx3,i1,i2,i3)
+
+ ! retrieve friction parameters
+ vo=structure(i3)%gammadot0
+ tauc=structure(i3)%stressexponent
+ friction=structure(i3)%friction
+ cohesion=structure(i3)%cohesion
- ! traction = sigma . n
- s=sig(i1,i2,i3)
- t=s .tdot. n
+ ! traction = sigma . n
+ s=sig(i1,i2,i3)
+ t=s .tdot. n
- ! signed normal component
- taun=SUM(t*n)
+ ! signed normal component
+ taun=SUM(t*n)
- ! absolute value of shear component
- ts=t-taun*n
- taus=SQRT(SUM(ts*ts))
+ ! absolute value of shear component
+ ts=t-taun*n
+ taus=SQRT(SUM(ts*ts))
- ! effective shear stress on fault plane
- tau=taus+friction*taun-cohesion
+ ! effective shear stress on fault plane
+ tau=MAX(0.d0,taus+friction*taun-cohesion)
- ! shear traction direction
- ts=ts/taus
+ ! rake direction test only if | rake | < 3*Pi
+ IF (SUM(ts*r).LT.0.d0 .AND. ABS(rake).LT.pi2*1.5d0) CYCLE
- ! 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
- ! creep rate
- slip=vo*2._8*sinh(tau/tauc)
+ ! creep rate
+ patch(j2,j3)%slip=vo*2._8*sinh(tau/tauc)
- ! strike-direction creep rate
- ss=slip*SUM(ts*sv)
+ ! shear traction direction
+ ts=ts/taus
- ! dip-direction creep rate
- ds=slip*SUM(ts*dv)
- END IF
- END IF
+ ! strike-direction creep rate
+ patch(j2,j3)%ss=patch(j2,j3)%slip*SUM(ts*sv)
- ! gather absolute and relative position, total,
- ! strike and dip slip in a single structure
- patch(j2,j3)=SLIPPATCH_STRUCT(x1,x2,x3,yr,zr,slip,ss,ds)
+ ! dip-direction creep rate
+ patch(j2,j3)%ds=patch(j2,j3)%slip*SUM(ts*dv)
END DO
END DO
diff -r 8d085e25c7e7 -r ed3436b7ac71 input.f90
--- a/input.f90 Wed Dec 14 15:04:59 2011 -0800
+++ b/input.f90 Tue Dec 27 16:04:21 2011 -0800
@@ -57,7 +57,7 @@ CONTAINS
!$ INTEGER :: omp_get_num_procs,omp_get_max_threads
REAL*8 :: dummy,dum1,dum2
REAL*8 :: minlength,minwidth
- TYPE(OPTION_S) :: opts(9)
+ TYPE(OPTION_S) :: opts(11)
INTEGER :: k,iostatus,i,e
@@ -69,15 +69,17 @@ CONTAINS
END IF
! parse the command line for options
- opts(1)=OPTION_S("no-proj-output",.FALSE.,CHAR(20))
- opts(2)=OPTION_S("no-relax-output",.FALSE.,CHAR(21))
- opts(3)=OPTION_S("no-txt-output",.FALSE.,CHAR(22))
- opts(4)=OPTION_S("no-vtk-output",.FALSE.,CHAR(23))
- opts(5)=OPTION_S("no-grd-output",.FALSE.,CHAR(24))
- opts(6)=OPTION_S("no-xyz-output",.FALSE.,CHAR(25))
- opts(7)=OPTION_S("no-stress-output",.FALSE.,CHAR(26))
- opts(8)=OPTION_S("dry-run",.FALSE.,CHAR(27))
- opts(9)=OPTION_S("help",.FALSE.,'h')
+ opts( 1)=OPTION_S("no-proj-output",.FALSE.,CHAR(20))
+ opts( 2)=OPTION_S("no-relax-output",.FALSE.,CHAR(21))
+ opts( 3)=OPTION_S("no-txt-output",.FALSE.,CHAR(22))
+ opts( 4)=OPTION_S("no-vtk-output",.FALSE.,CHAR(23))
+ opts( 5)=OPTION_S("no-grd-output",.FALSE.,CHAR(24))
+ opts( 6)=OPTION_S("no-xyz-output",.FALSE.,CHAR(25))
+ opts( 7)=OPTION_S("no-stress-output",.FALSE.,CHAR(26))
+ opts( 8)=OPTION_S("with-stress-output",.FALSE.,CHAR(27))
+ opts( 9)=OPTION_S("with-vtk-relax-output",.FALSE.,CHAR(28))
+ opts(10)=OPTION_S("dry-run",.FALSE.,CHAR(29))
+ opts(11)=OPTION_S("help",.FALSE.,'h')
DO
ch=getopt("h",opts)
@@ -107,6 +109,12 @@ CONTAINS
in%isoutputstress=.FALSE.
CASE(CHAR(27))
! option dry-run
+ in%isoutputstress=.TRUE.
+ CASE(CHAR(28))
+ ! option dry-run
+ in%isoutputvtkrelax=.TRUE.
+ CASE(CHAR(29))
+ ! option dry-run
in%isdryrun=.TRUE.
CASE('h')
! option help
@@ -129,16 +137,18 @@ CONTAINS
PRINT '(" [--no-vtk-output] [--no-xyz-output]")'
PRINT '("")'
PRINT '("options:")'
- PRINT '(" -h prints this message and aborts calculation")'
- PRINT '(" --dry-run abort calculation, only output geometry")'
- PRINT '(" --help prints this message and aborts calculation")'
- PRINT '(" --no-grd-output cancel output in GMT grd binary format")'
- PRINT '(" --no-proj-output cancel output in geographic projection")'
- PRINT '(" --no-relax-output cancel output of the postseismic contribution")'
- PRINT '(" --no-stress-output cancel output of stress tensor in any format")'
- PRINT '(" --no-txt-output cancel output in text format")'
- PRINT '(" --no-vtk-output cancel output in Paraview VTK format")'
- PRINT '(" --no-xyz-output cancel output in GMT xyz format")'
+ PRINT '(" -h prints this message and aborts calculation")'
+ PRINT '(" --dry-run abort calculation, only output geometry")'
+ PRINT '(" --help prints this message and aborts calculation")'
+ PRINT '(" --no-grd-output cancel output in GMT grd binary format")'
+ PRINT '(" --no-proj-output cancel output in geographic projection")'
+ PRINT '(" --no-relax-output cancel output of the postseismic contribution")'
+ PRINT '(" --no-stress-output cancel output of stress tensor in any format")'
+ PRINT '(" --no-txt-output cancel output in text format")'
+ PRINT '(" --no-vtk-output cancel output in Paraview VTK format")'
+ PRINT '(" --no-xyz-output cancel output in GMT xyz format")'
+ PRINT '(" --with-stress-output export stress tensor")'
+ PRINT '(" --with-vtk-relax-output export relaxation to VTK format")'
PRINT '("")'
PRINT '("description:")'
PRINT '(" Evaluates the deformation due to fault slip, surface loading")'
diff -r 8d085e25c7e7 -r ed3436b7ac71 makefile_imkl
--- a/makefile_imkl Wed Dec 14 15:04:59 2011 -0800
+++ b/makefile_imkl Tue Dec 27 16:04:21 2011 -0800
@@ -10,7 +10,7 @@
# to fit with your environment.
LIBPATH=-L/sw/lib -L/opt/intel/Compiler/11.1/084/Frameworks/mkl/lib/em64t/ -L/opt/intel/Compiler/11.1/084/lib/
-INCPATH=-I/sw/include -I/opt/intel/Compiler/11.1/084/Frameworks/mkl/include
+INCPATH=-I/sw/include -I/opt/intel/Compiler/11.1/084/Frameworks/mkl/include -Ibuild
LIBS=-lproj -lgmt -lnetcdf -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -openmp -lpthread
@@ -20,9 +20,9 @@ FC=ifort -openmp
FFLAGS=-cpp $(INCPATH) -zero -warn all
F77FLAGS=-zero
-CFLAGS=-I/sw/include
+CFLAGS=-I/sw/include -Ibuild
-OBJ = types.o mkl_dfti.o fourier.o green.o elastic3d.o friction3d.o viscoelastic3d.o writegrd4.2.o proj.o export.o getdata.o getopt.o input.o relax.o
+OBJ = types.o mkl_dfti.o fourier.o green.o elastic3d.o friction3d.o viscoelastic3d.o writegrd4.2.o writevtk.o proj.o export.o getdata.o getopt_m.o input.o relax.o
%.o : %.c
$(CC) $(CFLAGS) -c $^
diff -r 8d085e25c7e7 -r ed3436b7ac71 relax.f90
--- a/relax.f90 Wed Dec 14 15:04:59 2011 -0800
+++ b/relax.f90 Tue Dec 27 16:04:21 2011 -0800
@@ -217,6 +217,7 @@ PROGRAM relax
! arrays
REAL*4, DIMENSION(:,:,:), ALLOCATABLE :: v1,v2,v3,u1,u2,u3,gamma
+ REAL*4, DIMENSION(:,:,:), ALLOCATABLE :: u1r,u2r,u3r
REAL*4, DIMENSION(:,:), ALLOCATABLE :: t1,t2,t3
REAL*4, DIMENSION(:,:,:), ALLOCATABLE :: inter1,inter2,inter3
TYPE(TENSOR), DIMENSION(:,:,:), ALLOCATABLE :: tau,sig,moment
@@ -250,6 +251,16 @@ PROGRAM relax
t1(in%sx1+2,in%sx2),t2(in%sx1+2,in%sx2),t3(in%sx1+2,in%sx2), &
STAT=iostatus)
IF (iostatus>0) STOP "could not allocate memory"
+#ifdef VTK
+ IF (in%isoutputvtkrelax) THEN
+ ALLOCATE(u1r(in%sx1+2,in%sx2,in%sx3/2),u2r(in%sx1+2,in%sx2,in%sx3/2), &
+ u3r(in%sx1+2,in%sx2,in%sx3/2),STAT=iostatus)
+ IF (iostatus>0) STOP "could not allocate memory for VTK relax output"
+ u1r=0
+ u2r=0
+ u3r=0
+ END IF
+#endif
IF (in%isoutputrelax) THEN
ALLOCATE(inter1(in%sx1+2,in%sx2,2),inter2(in%sx1+2,in%sx2,2),inter3(in%sx1+2,in%sx2,2),STAT=iostatus)
@@ -383,6 +394,13 @@ PROGRAM relax
4,4,8,filename,title,name)
!CALL exportvtk_vectors_slice(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,8,8,filename)
END IF
+ IF (in%isoutputvtkrelax) THEN
+ filename=trim(in%wdir)//"/disp-relax-000.vtk"//char(0)
+ title="postseismic displacement vector field"//char(0)
+ name="displacement"//char(0)
+ CALL exportvtk_vectors_legacy(u1r,u2r,u3r,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3, &
+ 4,4,8,filename,title,name)
+ END IF
#endif
IF (ALLOCATED(in%ptsname)) THEN
CALL exportpoints(u1,u2,u3,sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
@@ -448,7 +466,7 @@ PROGRAM relax
CALL tensorfieldadd(moment,moment,in%sx1,in%sx2,in%sx3/2,c1=0._4,c2=0._4)
DO i=1,ITERATION_MAX
- IF (t > (in%interval+1e-6)) GOTO 100 ! proper exit
+ IF (t .GE. in%interval) GOTO 100 ! proper exit
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! predictor
@@ -573,8 +591,8 @@ PROGRAM relax
! nonlinear fault creep with rate-strengthening friction
IF (ALLOCATED(in%faultcreepstruc)) THEN
+
! use v1 as placeholders for the afterslip planes
- v1=0
DO k=1,in%np
! one may use optional arguments ...,VEL=v1) to convert
! fault slip to eigenstrain (scalar)
@@ -584,15 +602,13 @@ PROGRAM relax
sig,in%mu,in%faultcreepstruc,in%sx1,in%sx2,in%sx3/2, &
in%dx1,in%dx2,in%dx3,moment)
END DO
-
- ! update slip history
- CALL fieldadd(gamma,v1,in%sx1+2,in%sx2,in%sx3/2,c2=REAL(Dt))
! 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
! interseismic loading
@@ -644,6 +660,14 @@ PROGRAM relax
CALL sliceadd(inter2(:,:,1),v2,in%sx1+2,in%sx2,in%sx3,int(in%oz/in%dx3)+1,c2=REAL(Dt))
CALL sliceadd(inter3(:,:,1),v3,in%sx1+2,in%sx2,in%sx3,int(in%oz/in%dx3)+1,c2=REAL(Dt))
END IF
+
+#ifdef VTK
+ IF (in%isoutputvtkrelax) THEN
+ u1r=u1r+Dt*v1
+ u2r=u2r+Dt*v2
+ u3r=u3r+Dt*v3
+ END IF
+#endif
! time increment
t=t+Dt
@@ -749,6 +773,13 @@ PROGRAM relax
CALL exportvtk_vectors_legacy(v1,v2,v3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
8,8,16,filename,title,name)
!CALL exportvtk_vectors_slice(v1,v2,v3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,8,8,filename)
+ END IF
+ IF (in%isoutputvtkrelax) THEN
+ filename=trim(in%wdir)//"/disp-relax-"//digit//".vtk"//char(0)
+ title="postseismic displacement vector field"//char(0)
+ name="displacement"//char(0)
+ CALL exportvtk_vectors_legacy(u1r,u2r,u3r,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3, &
+ 4,4,8,filename,title,name)
END IF
#endif
diff -r 8d085e25c7e7 -r ed3436b7ac71 types.f90
--- a/types.f90 Wed Dec 14 15:04:59 2011 -0800
+++ b/types.f90 Tue Dec 27 16:04:21 2011 -0800
@@ -190,6 +190,7 @@ MODULE types
LOGICAL :: isoutputrelax=.TRUE.
LOGICAL :: isoutputtxt=.TRUE.
LOGICAL :: isoutputvtk=.TRUE.
+ LOGICAL :: isoutputvtkrelax=.TRUE.
LOGICAL :: isoutputgrd=.TRUE.
LOGICAL :: isoutputxyz=.TRUE.
LOGICAL :: isoutputstress=.TRUE.
diff -r 8d085e25c7e7 -r ed3436b7ac71 util/grdmap.sh
--- a/util/grdmap.sh Wed Dec 14 15:04:59 2011 -0800
+++ b/util/grdmap.sh Tue Dec 27 16:04:21 2011 -0800
@@ -27,7 +27,7 @@ done
done
# plotting vectors
-if [ -e "$U1" ]; then
+if [ -e $U1 ]; then
#echo $self": found "$U1": plotting vectors"
echo $self": Using VECTOR="$VECTOR", STEP="$ADX
@@ -60,7 +60,7 @@ fi
fi
#-Q0.20c/1.0c/0.4cn1.0c \
-psscale -O -Ba$PSSCALE/:$unit: -D2.0i/-0.8i/3.1i/0.2ih \
+psscale -O -B$PSSCALE/:$unit: -D2.0i/-0.8i/3.1i/0.2ih \
-C$colors $REVERT \
>> $PSFILE
@@ -116,7 +116,7 @@ do
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)}'`;;
+ 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;;
@@ -158,8 +158,8 @@ if [ "$tset" != "1" ]; then
if [ "$tset" != "1" ]; then
tick=`echo $bds | awk -F "/" '{s=1;print ($2-$1)/s/4}'`
fi
-HEIGHT=`echo $bds | awk -F "/" '{printf("%fi\n",($4-$3)/($2-$1)*4)}'`
if [ "$gset" != "1" ]; then
+ HEIGHT=`echo $bds | awk -F "/" '{printf("%fi\n",($4-$3)/($2-$1)*4)}'`
if [ "$rset" != "1" ]; then
PROJ="X4i/"$HEIGHT
else
@@ -167,7 +167,8 @@ if [ "$gset" != "1" ]; then
fi
AXIS=-Ba${tick}:"":/a${tick}:""::."$U3":WSne
else
- PROJ="M0/0/4i"
+ HEIGHT=4i
+ PROJ="M0/0/$HEIGHT"
AXIS=-B${tick}:"":/${tick}:""::."$U3":WSne
fi
@@ -208,7 +209,7 @@ if [ "$pset" != "1" ]; then
grd2cpt $U3 -C$cptfile -Z -T= > $colors
fi
else
- makecpt -C$cptfile -T$P -Z -D > $colors;
+ makecpt -C$cptfile -T$P -D > $colors;
fi
if [ "$oset" != "1" ]; then
ODIR=$WDIR
@@ -217,7 +218,7 @@ if [ -e $U1 ]; then
if [ "$vset" != "1" ]; then
MAX1=`grdinfo $U1 -C | awk '{if (sqrt($7^2)>sqrt($6^2)){print sqrt($7^2)}else{print sqrt($6^2)}}'`
MAX2=`grdinfo $U2 -C | awk '{if (sqrt($7^2)>sqrt($6^2)){print sqrt($7^2)}else{print sqrt($6^2)}}'`
- SIZE=`echo "$MAX1 $MAX2"| awk '{print ($1+$2)*0.95}'`
+ SIZE=`echo "$MAX1 $MAX2"| awk '{if (0==$1 && 0==$2){print 1}else{print ($1+$2)*0.95}}'`
VECTOR=$SIZE"c"
fi
if [ "$sset" != "1" ]; then
diff -r 8d085e25c7e7 -r ed3436b7ac71 writevtk.c
--- a/writevtk.c Wed Dec 14 15:04:59 2011 -0800
+++ b/writevtk.c Tue Dec 27 16:04:21 2011 -0800
@@ -206,12 +206,12 @@ void exportvtk_tensors_legacy_(sig,sx1,s
buffer[0]=(1u==endian)?sig[index+0]:swap(sig[index+0]);
buffer[1]=(1u==endian)?sig[index+1]:swap(sig[index+1]);
buffer[2]=(1u==endian)?sig[index+2]:swap(sig[index+2]);
- buffer[3]=(1u==endian)?sig[index+3]:swap(sig[index+3]);
- buffer[4]=(1u==endian)?sig[index+4]:swap(sig[index+4]);
- buffer[5]=(1u==endian)?sig[index+5]:swap(sig[index+5]);
- buffer[6]=(1u==endian)?sig[index+6]:swap(sig[index+6]);
- buffer[7]=(1u==endian)?sig[index+7]:swap(sig[index+7]);
- buffer[8]=(1u==endian)?sig[index+8]:swap(sig[index+8]);
+ buffer[3]=(1u==endian)?sig[index+1]:swap(sig[index+1]);
+ buffer[4]=(1u==endian)?sig[index+3]:swap(sig[index+3]);
+ buffer[5]=(1u==endian)?sig[index+4]:swap(sig[index+4]);
+ buffer[6]=(1u==endian)?sig[index+2]:swap(sig[index+2]);
+ buffer[7]=(1u==endian)?sig[index+4]:swap(sig[index+4]);
+ buffer[8]=(1u==endian)?sig[index+5]:swap(sig[index+5]);
// write buffer to disk
fwrite(buffer,36,1,funit);
More information about the CIG-COMMITS
mailing list