[cig-commits] commit: added the Savage and Prescott (1978) example
Mercurial
hg at geodynamics.org
Tue Sep 20 12:13:25 PDT 2011
changeset: 26:e44e659df0db
user: Sylvain Barbot <sylbar.vainbot at gmail.com>
date: Thu Aug 11 00:47:32 2011 -0700
files: examples/savage-prescott.sh export.f90 input.f90 relax.f90
description:
added the Savage and Prescott (1978) example
diff -r 19be26dce0d1 -r e44e659df0db examples/savage-prescott.sh
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/savage-prescott.sh Thu Aug 11 00:47:32 2011 -0700
@@ -0,0 +1,144 @@
+#!/bin/sh
+
+#
+# comments:
+#
+# linear viscoelastic relaxation following a strike-slip fault
+# followed by a second strike-slip event and a relaxation of
+# both stress perturbations.
+#
+# run this example with:
+#
+# ./run1.sh
+#
+# you can also abort computation to produce only geometry information with:
+#
+# ./run1.sh --dry-run
+#
+# viscous substrate starts at 3 times the seismogenic depth.
+# output at each computation step (Dt = -1)
+#
+# to visualize coseismic deformation (requires GRD output, or manual conversion to GRD format):
+#
+# map.sh -b -3/3/-3/3 output1/000
+#
+# to convert .xyz output in .grd output in post processing, type (requires GMT installed on your machine)
+#
+# xyz2grd.sh output1/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 output1/0{01,02,03,04,05,06,07,08,09,10}-relax
+#
+# type map.sh for a description of command-line options.
+# the command used to generate a map can be retrieve from the .ps file with
+#
+# tail -n 1 output1/000-plot.ps
+#
+# to visualize in 3-D with Paraview (www.paraview.org),
+# load the following files:
+#
+# output1/cgrid.vtp for the computational grid spatial extent
+# output1/rfaults-001.vtp for the slip distribution of the first event
+# output1/rfaults-002.vtp for the slip distribution of the second event
+# output1/disp-000.vtr for the coseismic displacement field
+# output1/linearlayer-001.vtp for the depth of brittle-ductile transition
+# output1/vel-001.vtr for the instantaneous velocity field
+#
+# note that we save the entire output to output1/in.param. this is convenient
+# to document how the simulation was computed.
+#
+# type relax --help to display options.
+#
+# modify the number of threads used with
+#
+# export OMP_NUM_THREADS=N with sh and ksh shells
+#
+# or
+#
+# setenv OMP_NUM_THREADS N with csh
+#
+
+count(){
+ echo $#
+}
+
+WDIR=./savage-prescott
+
+if [ ! -e $WDIR ]; then
+ echo adding directory $WDIR
+ mkdir $WDIR
+fi
+
+event="1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19"
+
+export OMP_NUM_THREADS=8
+
+time ../relax --no-proj-output --no-stress-output $* <<EOF | tee $WDIR/in.param
+# grid dimension (sx1,sx2,sx3)
+4 1024 512
+# sampling (dx1,dx2,dx3), smoothing (beta, nyquist)
+0.05 0.05 0.05 0.2 2
+# origin position (x0,y0) and rotation
+0 0 0
+# observation depth (displacement and stress)
+0 2
+# output directory
+$WDIR
+# lambda, mu, gamma (gamma = (1 - nu) rho g / mu)
+1 1 8.33e-4
+# time interval, (positive time step) or (negative skip, scaling)
+20 0.1
+# number of observation planes
+0
+# number of observation points
+0
+# number of stress observation segments
+0
+# number of prestress interfaces
+0
+# number of linear viscous interfaces
+1
+# n. depth gammadot0 cohesion
+ 1 4 1 0
+# number of linear weak zones
+0
+# number of nonlinear viscous interfaces
+0
+# number of fault creep interfaces
+0
+# number of inter-seismic strike-slip segments
+2
+# n. slip/time xs ys zs length width strike dip rake
+ 1 -1 -2 -7.0 0 4 3 0 90 0
+ 2 -1 -2 7.0 0 4 3 0 90 0
+# number of inter-seismic tensile segments
+0
+# number of events
+$(count $event)
+# number of coseismic strike-slip segments
+1
+# n. slip xs ys zs length width strike dip rake
+ 1 0 -6.4 0 0 12.8 1 0 90 0
+# number of coseismic tensile segments
+0
+# number of coseismic dilatation point sources
+0
+# number of surface loads
+0
+# generate a sequence of earthquakes
+`for e in $event; do
+ echo \# time of next coseismic event
+ echo $e
+ echo \# number of coseismic strike-slip segments
+ echo 1
+ echo \# n. slip xs ys zs length width strike dip rake
+ echo 1 1 -3.2 0 0 6.4 1 0 90 0
+ echo \# number of coseismic tensile segments
+ echo 0
+ echo \# number of coseismic dilatation point sources
+ echo 0
+ echo \# number of surface loads
+ echo 0
+done`
+EOF
diff -r 19be26dce0d1 -r e44e659df0db export.f90
--- a/export.f90 Wed Aug 10 19:04:03 2011 -0700
+++ b/export.f90 Thu Aug 11 00:47:32 2011 -0700
@@ -1326,6 +1326,8 @@ END SUBROUTINE exportcreep
REAL*8 :: friction
! traction components
REAL*8, DIMENSION(3) :: t,ts
+
+ IF (0.GE.nsop) RETURN
! double-quote character
q=char(34)
diff -r 19be26dce0d1 -r e44e659df0db input.f90
--- a/input.f90 Wed Aug 10 19:04:03 2011 -0700
+++ b/input.f90 Thu Aug 11 00:47:32 2011 -0700
@@ -769,7 +769,7 @@ CONTAINS
ALLOCATE(in%inter%s(in%inter%ns),in%inter%sc(in%inter%ns),STAT=iostatus)
IF (iostatus>0) STOP "could not allocate the source list"
PRINT 2000
- PRINT '(a)',"no. slip xs ys zs length width strike dip rake"
+ PRINT '(a)',"no. slip/time xs ys zs length width strike dip rake"
PRINT 2000
DO k=1,in%inter%ns
CALL getdata(iunit,dataline)
diff -r 19be26dce0d1 -r e44e659df0db relax.f90
--- a/relax.f90 Wed Aug 10 19:04:03 2011 -0700
+++ b/relax.f90 Thu Aug 11 00:47:32 2011 -0700
@@ -413,7 +413,7 @@ PROGRAM relax
#endif
#ifdef VTK
WRITE (digit4,'(I4.4)') oi-1
- IF (in%isoutputvtk) THEN
+ IF (in%isoutputvtk .AND. in%isoutputstress) THEN
filename=trim(in%wdir)//"/sigma-"//digit4//".vtk"
title="stress tensor field"
name="stress"
@@ -477,7 +477,7 @@ PROGRAM relax
END IF
#ifdef VTK
- IF (in%isoutputvtk) THEN
+ IF (in%isoutputvtk .AND. in%isoutputstress) THEN
WRITE (digit,'(I3.3)') oi-1
filename=trim(in%wdir)//"/power-"//digit//".vtk"
title="stress rate tensor field"
@@ -576,7 +576,7 @@ PROGRAM relax
IF ((in%inter%ns .GT. 0) .OR. (in%inter%nt .GT. 0)) THEN
! vectors v1,v2,v3 are not affected.
CALL dislocations(in%inter,in%lambda,in%mu,in%beta,in%sx1,in%sx2,in%sx3, &
- in%dx1,in%dx2,in%dx3,v1,v2,v3,t1,t2,t3,tau,factor=Dt,eigenstress=moment)
+ in%dx1,in%dx2,in%dx3,v1,v2,v3,t1,t2,t3,tau,eigenstress=moment)
END IF
v1=0;v2=0;v3=0;t1=0;t2=0;t3=0;
@@ -599,7 +599,7 @@ PROGRAM relax
#ifdef GRD_EQBF
IF (in%isoutputgrd) THEN
CALL exportgrd(v1,v2,v3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
- in%oz,in%x0,y0,in%wdir,oi,convention=3)
+ in%oz,in%x0,in%y0,in%wdir,oi,convention=3)
END IF
#endif
END IF
@@ -625,9 +625,10 @@ PROGRAM relax
IF (abs(t-in%events(e+1)%time) .LT. 1e-6) THEN
e=e+1
in%events(e)%i=i
+
PRINT '("coseismic event ",I3.3)', e
PRINT 0990
-
+
v1=0;v2=0;v3=0;t1=0;t2=0;t3=0;
CALL dislocations(in%events(e),in%lambda,in%mu, &
in%beta,in%sx1,in%sx2,in%sx3, &
@@ -792,7 +793,7 @@ CONTAINS
REAL*4, DIMENSION(sx1,sx2), INTENT(INOUT) :: t3
#endif
- INTEGER :: i1,i2,i3
+ INTEGER :: i,i1,i2,i3
DO i=1,e%nl
CALL shiftedindex(e%l(i)%x,e%l(i)%y,0._8,sx1,sx2,1,dx1,dx2,1._8,i1,i2,i3)
@@ -822,7 +823,7 @@ CONTAINS
TYPE(TENSOR), DIMENSION(:,:,:), INTENT(INOUT), OPTIONAL :: eigenstress
INTEGER :: i
- REAL*8 :: slip_factor=1._8
+ REAL*8 :: slip_factor
IF (PRESENT(factor)) THEN
slip_factor=factor
More information about the CIG-COMMITS
mailing list