[cig-commits] commit: export all ductile zones in a single .vtk file. Update the documentation to include lateral variations in viscoelasticity.
Mercurial
hg at geodynamics.org
Thu Oct 11 04:45:21 PDT 2012
changeset: 145:5b06b98308ff
tag: tip
user: Sylvain Barbot <sbarbot at caltech.edu>
date: Thu Oct 11 04:45:11 2012 -0700
files: latex/graphics/cover.pdf latex/relax.tex src/export.f90 src/input.f90
description:
export all ductile zones in a single .vtk file. Update the documentation to include lateral variations in viscoelasticity.
diff -r 703b8ccf773a -r 5b06b98308ff latex/graphics/cover.pdf
Binary file latex/graphics/cover.pdf has changed
diff -r 703b8ccf773a -r 5b06b98308ff latex/relax.tex
--- a/latex/relax.tex Mon Oct 08 10:45:18 2012 -0700
+++ b/latex/relax.tex Thu Oct 11 04:45:11 2012 -0700
@@ -42,14 +42,14 @@
\title{\bf Relax: Nonlinear Postseismic Relaxation\\in the Fourier Domain}
\author
-{User Manual\\
+{User's Manual\\
\\
-\normalsize{Sylvain Barbot (sbarbot at caltech.edu)}\\
-\normalsize{Division of Geological and Planetary Sciences, Caltech}\\
+\normalsize{Sylvain Barbot (sbarbot at ntu.edu.sg)}\\
+\normalsize{Earth Observatory of Singapore, Nanyang Technological University}\\
}
% Include the date command, but leave its argument blank.
-\date{Oct. 11th, 2011, last revision Dec. 28th, 2011}
+\date{Oct. 11th, 2011, last revision Oct. 11th, 2012}
\begin{document}
@@ -286,7 +286,7 @@ Many simulations imply the relaxation of
%
\begin{figure*}[h]
\centering\includegraphics[width=0.6\textwidth]{geologic_fault.pdf}
-\caption{Geologic model of slip on a fault. The fault is described by the position of its top tip, its length in the along-strike direction, width in the down-dip direction, its strike angle and its dip angle. The rake indicates the orientation of slip. For volumes (such as those used to define slabs and shear zones), an extra thickness parameter indicates how the rectangular volume is extruded from its center.}
+\caption{Geologic model of slip on a fault. The fault is described by the position of its top tip, its length in the along-strike direction, width in the down-dip direction, its strike angle and its dip angle. The rake indicates the orientation of slip.}
\label{fig:geologic_fault}
\end{figure*}
%
@@ -335,6 +335,38 @@ Let's look at how to define a viscosity
{\color{orange} 5 70.0 0.2 0.0}
\end{alltt}
Here, the lower crust relaxes with a Maxwell time of $t_m=10\,$yr, and the asthenosphere with a Maxwell time of $t_m=5\,$yr. There is no need to specify the viscosity in the upper 15\,km, which is assumed infinite.
+
+
+
+
+\subsection{Lateral variations of viscosity}
+Subducting slabs and shear zones below a transform boundary are examples of major structural elements that introduce lateral variations in viscosity in the lithosphere. Including these elements in the Earth structure is done using ductile zones, which are rectangular regions where the viscosity is modified from the background, depth-dependent, value. The volume where the viscosity changes is defined similarly to a fault (using start position, length, width, strike and dip), but a additional thickness parameter indicates how the volume is extruded from the central plane (Fig.~\ref{fig:geologic_zone}).
+
+%
+\begin{figure*}[h]
+\centering\includegraphics[width=0.6\textwidth]{geologic_zone.pdf}
+\caption{Geometry of a ductile zone anomaly. The ductile zone is described by the position of the tip of the central plane (length in the strike direction, width in the dip direction, strike and dip angles). The thickness parameter indicates how the rectangular volume is extruded from its center plane. The $\Delta\dot{\gamma}_0$ parameter adds to the background value of the fluidity.}
+\label{fig:geologic_zone}
+\end{figure*}
+%
+
+
+Let's look at how to define a viscosity profile where the depth of the brittle-ductile transition changes across a geological boundary. We start with an elastic plate overriding a viscoelastic half space with a uniform brittle-ductile transition depth of 15\,km with a Maxwell time of $t_m=1\,$yr.
+\begin{alltt}
+...
+# number of linear viscous interfaces
+{\color{orange}1}
+# no depth gammadot0 cohesion
+{\color{orange} 1 15.0 1 0.0}
+# number of linear ductile shear zones
+{\color{orange}1}
+# no dgammadot0 x1 x2 x3 length width thickness strike dip
+{\color{orange} 1 -1 0 0 10 200 5 200 90 90}
+...
+\end{alltt}
+We then add a large volume abutting the center of the computational volume, where we modify the fluidity from $\dot{\gamma}_0=1\,$yr$^{-1}$ to $\dot{\gamma}_0=0$, using $\Delta\dot{\gamma}_0=-1\,$yr$^{-1}$. As a null fluidity corresponds to elastic behavior, the brittle-ductile transition occurs 5\,km deeper going west to east across the $x_1$ (north-south) axis. The geometry of a slab can be obtained by modifying the strike and dip angles. Changing the strike of the ductile zone from 90 to 0$^\circ$ would make the transition occur from south to north. If the background and anomalous fluidities add to a negative, unphysical, value, the sum is corrected to null, and the result is an effective elastic property. Realistic structures can be accounted for from structural data using a large quantity of ductile anomalies tuned to observations.
+
+
\pagebreak
\section{Examples}
@@ -543,7 +575,7 @@ 0
0
...
\end{alltt}
-Note that the time step asked by the program is only the output time step, that is, at what time to output the result. Internally, an optimal time step is always evaluated and the real time step used for the calculation is the smallest of either the internally evaluated value or the time step required to produce an output every multiple of the required time step. The internally evaluated value is always $\Delta t^{\text{internal}}=t_m/10$. Prescribing the output time step can only reduce the computational time step. Note that is a background viscosity has been established, the program requests so called ``ductile zones", which are finite volumes where the background viscosity is perturbed. These volumes are described similarly as for a fault, but with an extra ``thickness" parameter indicating how the volumes extrudes from the center plane. The standard output of this simulation looks as follows
+Note that the time step asked by the program is only the output time step, that is, at what time to output the result. Internally, an optimal time step is always evaluated and the real time step used for the calculation is the smallest of either the internally evaluated value or the time step required to produce an output every multiple of the required time step. The internally evaluated value is always $\Delta t^{\text{internal}}=t_m/10$. Prescribing the output time step can only reduce the computational time step. Note that is a background viscosity has been established, the program requests so called ``ductile zones", which are finite volumes where the background viscosity is perturbed. These volumes are described similarly to for a fault, but with an extra ``thickness" parameter indicating how the volumes extrudes from the center plane. The standard output of this simulation looks as follows
\begin{alltt}
coseismic event 001
I | Dt | tm(ve) | tm(pl) | tm(as) | t/tmax | power | C:E^i |
diff -r 703b8ccf773a -r 5b06b98308ff src/export.f90
--- a/src/export.f90 Mon Oct 08 10:45:18 2012 -0700
+++ b/src/export.f90 Thu Oct 11 04:45:11 2012 -0700
@@ -2451,6 +2451,154 @@ END SUBROUTINE exportcreep_vtk
END SUBROUTINE exportvtk_brick
!------------------------------------------------------------------
+ !> subroutine ExportVTK_AllBricks
+ !! creates a .vtp file (in the VTK PolyData XML format) containing
+ !! a brick (3d rectangle, cuboid).
+ !!
+ !! \author sylvain barbot 06/24/09 - original form
+ !------------------------------------------------------------------
+ SUBROUTINE exportvtk_allbricks(swz,weakzones,filename)
+ INTEGER, INTENT(IN) :: swz
+ TYPE(WEAK_STRUCT), DIMENSION(swz), INTENT(IN) :: weakzones
+ CHARACTER(256), INTENT(IN) :: filename
+
+ INTEGER :: iostatus,i
+ CHARACTER :: q
+
+ REAL*8 :: x1,x2,x3,L,W,T,strike,dip
+ REAL*8 :: cstrike,sstrike,cdip,sdip
+ REAL*8, DIMENSION(3) :: s,d,n
+
+ ! double-quote character
+ q=char(34)
+
+ OPEN (UNIT=15,FILE=filename,IOSTAT=iostatus,FORM="FORMATTED")
+ IF (iostatus>0) THEN
+ WRITE_DEBUG_INFO
+ PRINT '(a)', filename
+ STOP "could not open file for export in ExportVTK_Brick"
+ END IF
+
+ WRITE (15,'("<?xml version=",a,"1.0",a,"?>")') q,q
+ WRITE (15,'("<VTKFile type=",a,"PolyData",a," version=",a,"0.1",a,">")') q,q,q,q
+ WRITE (15,'(" <PolyData>")')
+
+ DO i=1,swz
+ x1=weakzones(i)%x
+ x2=weakzones(i)%y
+ x3=weakzones(i)%z
+ L=weakzones(i)%length
+ W=weakzones(i)%width
+ T=weakzones(i)%thickness
+ strike=weakzones(i)%strike
+ dip=weakzones(i)%dip
+
+ 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
+
+ ! surface normal vector components
+ n(1)=+cdip*cstrike
+ n(2)=-cdip*sstrike
+ n(3)=-sdip
+
+ WRITE (15,'(" <Piece NumberOfPoints=",a,"8",a," NumberOfPolys=",a,"1",a,">")'),q,q,q,q
+ WRITE (15,'(" <Points>")')
+ WRITE (15,'(" <DataArray type=",a,"Float32",a, &
+ & " Name=",a,"Weak Zone",a, &
+ & " NumberOfComponents=",a,"3",a, &
+ & " format=",a,"ascii",a,">")'),q,q,q,q,q,q,q,q
+
+ ! fault edge coordinates
+ WRITE (15,'(24ES11.2)') &
+ x1-d(1)*W/2-s(1)*L/2-n(1)*T/2.d0, x2-d(2)*W/2-s(2)*L/2-n(2)*T/2.d0, x3-d(3)*W/2-s(3)*L/2-n(3)*T/2.d0, &
+ x1-d(1)*W/2+s(1)*L/2-n(1)*T/2.d0, x2-d(2)*W/2+s(2)*L/2-n(2)*T/2.d0, x3-d(3)*W/2+s(3)*L/2-n(3)*T/2.d0, &
+ x1+d(1)*W/2+s(1)*L/2-n(1)*T/2.d0, x2+d(2)*W/2+s(2)*L/2-n(2)*T/2.d0, x3+d(3)*W/2+s(3)*L/2-n(3)*T/2.d0, &
+ x1+d(1)*W/2-s(1)*L/2-n(1)*T/2.d0, x2+d(2)*W/2-s(2)*L/2-n(2)*T/2.d0, x3+d(3)*W/2-s(3)*L/2-n(3)*T/2.d0, &
+ x1+d(1)*W/2-s(1)*L/2+n(1)*T/2.d0, x2+d(2)*W/2-s(2)*L/2+n(2)*T/2.d0, x3+d(3)*W/2-s(3)*L/2+n(3)*T/2.d0, &
+ x1-d(1)*W/2-s(1)*L/2+n(1)*T/2.d0, x2-d(2)*W/2-s(2)*L/2+n(2)*T/2.d0, x3-d(3)*W/2-s(3)*L/2+n(3)*T/2.d0, &
+ x1-d(1)*W/2+s(1)*L/2+n(1)*T/2.d0, x2-d(2)*W/2+s(2)*L/2+n(2)*T/2.d0, x3-d(3)*W/2+s(3)*L/2+n(3)*T/2.d0, &
+ x1+d(1)*W/2+s(1)*L/2+n(1)*T/2.d0, x2+d(2)*W/2+s(2)*L/2+n(2)*T/2.d0, x3+d(3)*W/2+s(3)*L/2+n(3)*T/2.d0
+
+ WRITE (15,'(" </DataArray>")')
+ WRITE (15,'(" </Points>")')
+ WRITE (15,'(" <Polys>")')
+ WRITE (15,'(" <DataArray type=",a,"Int32",a, &
+ & " Name=",a,"connectivity",a, &
+ & " format=",a,"ascii",a, &
+ & " RangeMin=",a,"0",a, &
+ & " RangeMax=",a,"6",a,">")'), q,q,q,q,q,q,q,q,q,q
+ WRITE (15,'("7 4 5 6 7 4 3 2 7 2 1 6")')
+ WRITE (15,'(" </DataArray>")')
+ WRITE (15,'(" <DataArray type=",a,"Int32",a, &
+ & " Name=",a,"offsets",a, &
+ & " format=",a,"ascii",a, &
+ & " RangeMin=",a,"12",a, &
+ & " RangeMax=",a,"12",a,">")'), q,q,q,q,q,q,q,q,q,q
+ WRITE (15,'(" 12")')
+ WRITE (15,'(" </DataArray>")')
+ WRITE (15,'(" </Polys>")')
+ WRITE (15,'(" </Piece>")')
+
+ WRITE (15,'(" <Piece NumberOfPoints=",a,"8",a," NumberOfPolys=",a,"1",a,">")'),q,q,q,q
+ WRITE (15,'(" <Points>")')
+ WRITE (15,'(" <DataArray type=",a,"Float32",a, &
+ & " Name=",a,"Weak Zone",a, &
+ & " NumberOfComponents=",a,"3",a, &
+ & " format=",a,"ascii",a,">")'),q,q,q,q,q,q,q,q
+
+ ! fault edge coordinates
+ WRITE (15,'(24ES11.2)') &
+ x1-d(1)*W/2.d0-s(1)*L/2.d0-n(1)*T/2.d0, x2-d(2)*W/2.d0-s(2)*L/2.d0-n(2)*T/2.d0,x3-d(3)*W/2-s(3)*L/2-n(3)*T/2.d0, &
+ x1-d(1)*W/2.d0+s(1)*L/2.d0-n(1)*T/2.d0, x2-d(2)*W/2.d0+s(2)*L/2.d0-n(2)*T/2.d0,x3-d(3)*W/2+s(3)*L/2-n(3)*T/2.d0, &
+ x1+d(1)*W/2.d0+s(1)*L/2.d0-n(1)*T/2.d0, x2+d(2)*W/2.d0+s(2)*L/2.d0-n(2)*T/2.d0,x3+d(3)*W/2+s(3)*L/2-n(3)*T/2.d0, &
+ x1+d(1)*W/2.d0-s(1)*L/2.d0-n(1)*T/2.d0, x2+d(2)*W/2.d0-s(2)*L/2.d0-n(2)*T/2.d0,x3+d(3)*W/2-s(3)*L/2-n(3)*T/2.d0, &
+ x1+d(1)*W/2.d0-s(1)*L/2.d0+n(1)*T/2.d0, x2+d(2)*W/2.d0-s(2)*L/2.d0+n(2)*T/2.d0,x3+d(3)*W/2-s(3)*L/2+n(3)*T/2.d0, &
+ x1-d(1)*W/2.d0-s(1)*L/2.d0+n(1)*T/2.d0, x2-d(2)*W/2.d0-s(2)*L/2.d0+n(2)*T/2.d0,x3-d(3)*W/2-s(3)*L/2+n(3)*T/2.d0, &
+ x1-d(1)*W/2.d0+s(1)*L/2.d0+n(1)*T/2.d0, x2-d(2)*W/2.d0+s(2)*L/2.d0+n(2)*T/2.d0,x3-d(3)*W/2+s(3)*L/2+n(3)*T/2.d0, &
+ x1+d(1)*W/2.d0+s(1)*L/2.d0+n(1)*T/2.d0, x2+d(2)*W/2.d0+s(2)*L/2.d0+n(2)*T/2.d0,x3+d(3)*W/2+s(3)*L/2+n(3)*T/2.d0
+
+ WRITE (15,'(" </DataArray>")')
+ WRITE (15,'(" </Points>")')
+ WRITE (15,'(" <Polys>")')
+ WRITE (15,'(" <DataArray type=",a,"Int32",a, &
+ & " Name=",a,"connectivity",a, &
+ & " format=",a,"ascii",a, &
+ & " RangeMin=",a,"0",a, &
+ & " RangeMax=",a,"7",a,">")'), q,q,q,q,q,q,q,q,q,q
+ WRITE (15,'("0 1 2 3 0 5 4 3 0 1 6 5")')
+ WRITE (15,'(" </DataArray>")')
+ WRITE (15,'(" <DataArray type=",a,"Int32",a, &
+ & " Name=",a,"offsets",a, &
+ & " format=",a,"ascii",a, &
+ & " RangeMin=",a,"12",a, &
+ & " RangeMax=",a,"12",a,">")'), q,q,q,q,q,q,q,q,q,q
+ WRITE (15,'(" 12")')
+ WRITE (15,'(" </DataArray>")')
+ WRITE (15,'(" </Polys>")')
+ WRITE (15,'(" </Piece>")')
+
+ END DO
+
+ WRITE (15,'(" </PolyData>")')
+ WRITE (15,'("</VTKFile>")')
+
+ CLOSE(15)
+
+ END SUBROUTINE exportvtk_allbricks
+
+ !------------------------------------------------------------------
!> subroutine ExportVTK_Vectors
!! creates a .vtr file (in the VTK Rectilinear XML format)
!! containing a vector field.
diff -r 703b8ccf773a -r 5b06b98308ff src/input.f90
--- a/src/input.f90 Mon Oct 08 10:45:18 2012 -0700
+++ b/src/input.f90 Thu Oct 11 04:45:11 2012 -0700
@@ -567,21 +567,26 @@ CONTAINS
in%linearweakzone(k)%strike,in%linearweakzone(k)%dip, &
dummy,in%x0,in%y0,in%rot)
- WRITE (digit,'(I3.3)') k
+ WRITE (digit,'(I3.3)') k
#ifdef VTK
- ! export the ductile zone in VTK format
- rffilename=trim(in%wdir)//"/weakzone-"//digit//".vtp"
- CALL exportvtk_brick(in%linearweakzone(k)%x,in%linearweakzone(k)%y,in%linearweakzone(k)%z, &
- in%linearweakzone(k)%length,in%linearweakzone(k)%width,in%linearweakzone(k)%thickness, &
- in%linearweakzone(k)%strike,in%linearweakzone(k)%dip,rffilename)
+ ! export the ductile zone in VTK format
+ rffilename=trim(in%wdir)//"/weakzone-"//digit//".vtp"
+ CALL exportvtk_brick(in%linearweakzone(k)%x,in%linearweakzone(k)%y,in%linearweakzone(k)%z, &
+ in%linearweakzone(k)%length,in%linearweakzone(k)%width,in%linearweakzone(k)%thickness, &
+ in%linearweakzone(k)%strike,in%linearweakzone(k)%dip,rffilename)
#endif
- ! export the ductile zone in GMT .xy format
- rffilename=trim(in%wdir)//"/weakzone-"//digit//".xy"
- CALL exportxy_brick(in%linearweakzone(k)%x,in%linearweakzone(k)%y,in%linearweakzone(k)%z, &
- in%linearweakzone(k)%length,in%linearweakzone(k)%width,in%linearweakzone(k)%thickness, &
- in%linearweakzone(k)%strike,in%linearweakzone(k)%dip,rffilename)
+ ! export the ductile zone in GMT .xy format
+ rffilename=trim(in%wdir)//"/weakzone-"//digit//".xy"
+ CALL exportxy_brick(in%linearweakzone(k)%x,in%linearweakzone(k)%y,in%linearweakzone(k)%z, &
+ in%linearweakzone(k)%length,in%linearweakzone(k)%width,in%linearweakzone(k)%thickness, &
+ in%linearweakzone(k)%strike,in%linearweakzone(k)%dip,rffilename)
END DO
+#ifdef VTK
+ ! export the ductile zone in VTK format
+ rffilename=trim(in%wdir)//"/weakzones-linear.vtp"
+ CALL exportvtk_allbricks(in%nlwz,in%linearweakzone,rffilename)
+#endif
END IF
END IF ! end linear viscous
@@ -706,6 +711,11 @@ CONTAINS
in%nonlinearweakzone(k)%strike, &
in%nonlinearweakzone(k)%dip,rffilename)
END DO
+#ifdef VTK
+ ! export the ductile zone in VTK format
+ rffilename=trim(in%wdir)//"/weakzones-nonlinear.vtp"
+ CALL exportvtk_allbricks(in%nnlwz,in%nonlinearweakzone,rffilename)
+#endif
END IF
END IF ! end nonlinear viscous
More information about the CIG-COMMITS
mailing list