[cig-commits] commit: More edits for reviewer comments. Started work on response to JGR reviews.
Mercurial
hg at geodynamics.org
Tue Feb 5 14:21:26 PST 2013
changeset: 160:ea734d8e36ef
tag: tip
user: Brad Aagaard <baagaard at usgs.gov>
date: Tue Feb 05 14:21:20 2013 -0800
files: faultRup.tex response_jgr.tex response_usgs.tex
description:
More edits for reviewer comments. Started work on response to JGR reviews.
diff -r 4c38d2c89757 -r ea734d8e36ef faultRup.tex
--- a/faultRup.tex Wed Jan 23 07:41:39 2013 -0800
+++ b/faultRup.tex Tue Feb 05 14:21:20 2013 -0800
@@ -123,9 +123,9 @@ thousands of years of strain accumulatio
thousands of years of strain accumulation between earthquakes. The
complexity of the many physical processes operating over this vast
range of scales leads most researchers to focus on a narrow space-time
-window in order to isolate just one or a few processes; the limited
-spatial and temporal coverage of observations also often justifies
-this narrow focus.
+window to isolate just one or a few processes; the limited spatial and
+temporal coverage of observations also often justifies this narrow
+focus.
Researchers have recognized for some time, though, that interseismic
deformation and fault interactions influence earthquake rupture
@@ -137,21 +137,20 @@ the static coseismic slip
the static coseismic slip
\citep{Reilinger:etal:2000,Pollitz:etal:2001,Langbein:etal:2006,Chlieh:etal:2007}.
Likewise, studies of rapid deformation associated with earthquake
-rupture propagation often approximate the loading of the crust at the beginning of a
-rupture
+rupture propagation often approximate the loading of the crust at the
+beginning of a rupture
\citep{Mikumo:etal:1998,Harris:Day:1999,Aagaard:etal:BSSA:2001,Peyrat:etal:2001,Oglesby:Day:2001,Dunham:Archuleta:2004}. Numerical
seismicity models that attempt to model multiple earthquake cycles,
generally simplify not only the fault loading and rupture propagation
-but also the physical properties in order to make the calculations
-tractable
+but also the physical properties to make the calculations tractable
\citep{Ward:1992,Robinson:Benites:1995,Hillers:etal:2006,Rundle:etal:2006,Pollitz:Schwartz:2008,Dieterich:Richards-Dinger:2010}.
Some dynamic spontaneous rupture modeling studies have
-attempted to examine a broader space-time window in order to remove
+attempted to examine a broader space-time window to remove
simplifying assumptions and more accurately capture the complex
interactions over the earthquake cycle. For example,
\citet{Duan:Oglesby:2005} simulated multiple earthquake cycles on a
-fault with a bend in order to capture the spatial variation in the
+fault with a bend to capture the spatial variation in the
stress field around the bend, which they found to have a strong role
in determining whether a rupture would propagate past the bend. By
spinning up the model over many earthquake cycles, they obtained a
@@ -181,7 +180,7 @@ the amount of aseismic slip.
the amount of aseismic slip.
Collectively, these studies suggest a set of desirable features for
-models of the earthquake cycle in order to capture both the slow
+models of the earthquake cycle to capture both the slow
deformation associated with interseismic behavior and the rapid
deformation associated with earthquake rupture propagation. These
features include the general capabilities of modeling elasticity with
@@ -313,7 +312,7 @@ codes
codes
\citep{Andrews:1999,Bizzarri:Cocco:2005,Day:etal:2005,Duan:Oglesby:2005,Dalguer:Day:2007,Moczo:etal:2007}
(ADD CITATION TO MELOSH AND RAEFSKY BSSA 1981 AS WELL, BUT CALL IT SPLIT NODES?),
-but differs from imposing fault slip via double-couple point
+but differs from imposing fault slip via double couple point
sources. The domain decomposition approach treats the fault surface as
a frictional contact interface, and the tractions correspond directly
to the Lagrange multipliers needed to satisfy the constraint equation
@@ -544,14 +543,15 @@ basis functions. This makes the traditio
basis functions. This makes the traditional LBB stability criterion
(:TODO: ADD FEM REFERENCE) trivial to satisfy by choosing the space of
Lagrange multipliers to be exactly the space of displacements,
-restricted to the fault. In simple terms in order to specify the
-problem we need to know the distance between any pair of fault
-vertices $(v^+,v^-)$, which can be expressed as a displacement.
+restricted to the fault. In simple terms to specify the problem we
+need to know the distance between any pair of vertices spanning the
+fault, which can be expressed as a relative displacement (i.e., fault
+slip).
% ------------------------------------------------------------------
\subsection{Dynamic Simulations}
-In dynamic simulations we include the inertial term in order to
+In dynamic simulations we include the inertial term to
resolve the propagation of seismic waves, with an intended focus on
applications to earthquake physics and ground-motion simulations. The
general form of the system Jacobian remains the same as in
@@ -821,7 +821,7 @@ finite-element operations.
In our domain decomposition approach, the finite-element mesh includes
the fault as an interior surface. This forces alignment of the element
-faces along the fault. In order to impose a given fault slip, as in
+faces along the fault. To impose a given fault slip as in
equation~(\ref{eqn:fault:disp}), we must represent the displacement on
both sides of the fault for any vertex on the fault. One option is to
designate ``fault vertices'' which possess twice as many displacement
@@ -920,7 +920,7 @@ matrix-vector multiplication is scalable
matrix-vector multiplication is scalable via parallel processing, this
is the method of choice for parallel simulation. However, for most
physically relevant problems, the Krylov solver requires a
-preconditioner in order to accelerate convergence. While generic
+preconditioner to accelerate convergence. While generic
preconditioners exist~\citep{Saad03,Smith:etal:1996}, the method must
often be specialized to a particular problem. In this section we
describe a preconditioner specialized to our formulation for fault
@@ -974,16 +974,16 @@ The elastic submatrix $\bm{K}$, in the a
The elastic submatrix $\bm{K}$, in the absence of boundary conditions,
has three translational and three rotational null modes. These are
provided to the algebraic multigrid (AMG) preconditioner, such as the
-ML library \citep{ML:users:guide} or the PETSc GAMG preconditioner, in
-order to assure an accurate coarse grid solution. AMG mimics the
-action of traditional geometric multigrid, but it generates coarse
-level operators and interpolation matrices using only the system
-matrix, treated as a weighted graph, rather than a separate
-description of the problem geometry, such as a mesh. We split the
-elastic block from the fault block and also manage the Schur
-complements. In this way, all block preconditioners, including those
-nested with multigrid, can be controlled from the options file without
-recompilation or special code.
+ML library \citep{ML:users:guide} or the PETSc GAMG preconditioner, to
+assure an accurate coarse grid solution. AMG mimics the action of
+traditional geometric multigrid, but it generates coarse level
+operators and interpolation matrices using only the system matrix,
+treated as a weighted graph, rather than a separate description of the
+problem geometry, such as a mesh. We split the elastic block from the
+fault block and also manage the Schur complements. In this way, all
+block preconditioners, including those nested with multigrid, can be
+controlled from the options file without recompilation or special
+code.
We now turn our attention to evaluating the fault portion of the
preconditioning matrix associated with the Lagrange multipliers, since
@@ -1029,7 +1029,7 @@ In dynamic simulations the Courant-Fride
\citep{Courant:etal:1967} controls the stability of the explicit time
integration. In most dynamic problems this dictates a relatively small
time step so that a typical simulation involves tens of thousands of
-time steps. Hence, we want a very efficient solver in order to run
+time steps. Hence, we want a very efficient solver to run
dynamic simulations in a reasonable amount of time.
% Lumped Jacobian
@@ -1048,7 +1048,7 @@ reduced by storing the diagonal of the m
reduced by storing the diagonal of the matrix as a vector rather
than as a sparse matrix. However, the block structure of our Jacobian
matrix, with the fault slip constraints occupying off-diagonal blocks,
-requires a two step approach in order to solve the linear system
+requires a two step approach to solve the linear system
of equations without forming a sparse matrix.
First, we eliminate the off-diagonal entries in each block of the
@@ -1126,7 +1126,7 @@ taking advantage of the fact that we con
taking advantage of the fact that we construct $\bm{K}$ so that it
is diagonal.
-We next compute the increment in the Lagrange multipliers in order to
+We next compute the increment in the Lagrange multipliers to
correct this initial solution so that the true residual is zero.
Making use of the initial residual, the expression for the true
residual is
@@ -1229,11 +1229,11 @@ the problem size (number of DOF) for the
the problem size (number of DOF) for the two different cell types
(hexahedra and tetrahedra) are nearly the same (within 2\%). The suite
of simulations examine increasing larger problem sizes as we increase
-the number of processes, with $7.8\times 10^4$ DOF for 1 process up to
-$7.1\times 10^6$ DOF for 96 processes. The corresponding
-discretization sizes are 2033 m to 437 m for the hexahedral meshes and
-2326 m to 712 m for the tetrahedral meshes.
-Figure~\ref{fig:solvertest:mesh} shows the 1846 m resolution
+the number of processes (with one process per core), with $7.8\times
+10^4$ DOF for 1 process up to $7.1\times 10^6$ DOF for 96
+processes. The corresponding discretization sizes are 2033 m to 437 m
+for the hexahedral meshes and 2326 m to 712 m for the tetrahedral
+meshes. Figure~\ref{fig:solvertest:mesh} shows the 1846 m resolution
tetrahedral mesh. As we will see in
Section~\ref{sec:verification:quasi-static}, the hexahedral mesh for a
given resolution in a quasi-static problem is slightly more accurate,
@@ -1298,9 +1298,9 @@ weighting the cohesive cells the same as
weighting the cohesive cells the same as conventional bulk cells while
partitioning. In this performance benchmark matrix-vector
multiplication (the PETSc \texttt{MatMult} function) has a load
-imbalance of up to 20\% on 96 processors. The cell partition balances
+imbalance of up to 20\% on 96 cores The cell partition balances
the number of cells across the processes using ParMetis
-\citep{Karypis:etal:1999} in order to achieve good balance for the
+\citep{Karypis:etal:1999} to achieve good balance for the
finite element integration. This does not take into account a
reduction in the number of DOF associated with constraints from
Dirichlet boundary conditions or the additional DOF associated with
@@ -1312,19 +1312,20 @@ We evaluate the parallel performance via
We evaluate the parallel performance via a weak scaling
criterion. That is, we run simulations on various numbers of
processors/cores with an increase in the problem size as the number of
-processes increases to maintain the same workload (e.g., number of
-cells and number of DOF) for each core. In ideal weak scaling the time
-for the various stages of the simulation is independent of the number
-of processes. For this performance benchmark we use the entire suite
-of hexahedral and tetrahedral meshes described earlier that range in
-size from $7.8\times 10^4$ DOF (1 process) to $7.1\times 10^6$ DOF (96
-processes). We employ the AMG preconditioner for the elasticity
-submatrix and our custom preconditioner for the Lagrange multipliers
-submatrix. We ran the simulations on Lonestar at the Texas Advanced
-Computing Center. Lonestar is comprised of 1888 compute nodes
-connected by QDR Infiniband in a fat-tree topology, where each compute
-node consisted of two six-core Intel Xeon E5650 processors with 24 GB
-of RAM. Simulations run on twelve or fewer cores were run on a single
+processes (with one process per core) increases to maintain the same
+workload (e.g., number of cells and number of DOF) for each core. In
+ideal weak scaling the time for the various stages of the simulation
+is independent of the number of processes. For this performance
+benchmark we use the entire suite of hexahedral and tetrahedral meshes
+described earlier that range in size from $7.8\times 10^4$ DOF (1
+process) to $7.1\times 10^6$ DOF (96 processes). We employ the AMG
+preconditioner for the elasticity submatrix and our custom
+preconditioner for the Lagrange multipliers submatrix. We ran the
+simulations on Lonestar at the Texas Advanced Computing
+Center. Lonestar is comprised of 1888 compute nodes connected by QDR
+Infiniband in a fat-tree topology, where each compute node consisted
+of two six-core Intel Xeon E5650 processors with 24 GB of
+RAM. Simulations run on twelve or fewer cores were run on a single
compute node with processes distributed across processors and then
cores. For example, the two process simulation used one core on each
of two processors. In addition to algorithm bottlenecks, runtime
@@ -1338,9 +1339,9 @@ given by the \texttt{VecAXPY}, \texttt{V
given by the \texttt{VecAXPY}, \texttt{VecMAXPY} and \texttt{VecMDot}
operations reported in the log summary \citep{PETSc:manual}. These
operations are limited by available memory bandwidth rather than the
-rate at which a processor can perform floating points operations. From
+rate at which a processor or core can perform floating points operations. From
Table~\ref{tab:solvertest:memory:events}, we see that we saturate the
-memory bandwidth using two processes per processor, since scaling
+memory bandwidth using two processes (cores) per processor, since scaling
plateaus from 2 to 4 processes, but shows good scaling from 12 to 24
processes. This lack of memory bandwidth will depress overall
performance, but should not affect the inter-node scaling of the
@@ -1391,15 +1392,19 @@ and reverse viscoelastic simulations and
and reverse viscoelastic simulations and various exercises in the
suite of dynamic spontaneous rupture benchmarks developed by the
Southern California Earthquake Center (SCEC) and the United States
-Geological Survey \citep{Harris:etal:SRL:2009}. In this section we
-focus on two benchmarks that test different scientific applications:
-quasi-static relaxation of a Maxwell viscoelastic material subjected
-to multiple earthquake cycles involving slip and steady creep on a
-vertical strike-slip fault \citep{Savage:Prescott:1978} and supershear
-dynamic spontaneous rupture of a 60 degree dipping normal fault in a
-Drucker-Prager elastoplastic medium. This second benchmark corresponds
-to benchmark TPV13 in the SCEC suite of dynamic spontaneous rupture
-benchmarks \citep{Harris:etal:SRL:2011}.
+Geological Survey \citep{Harris:etal:SRL:2009}. The mesh generation
+and simulation parameter files for many of the benchmarks, including
+those discussed here, are available from the CIG subversion repository
+(http://geodynamics.org/svn/cig/short/3D/PyLith/benchmarks/trunk/). In
+this section we focus on two benchmarks that test different scientific
+applications: quasi-static relaxation of a Maxwell viscoelastic
+material subjected to multiple earthquake cycles involving slip and
+steady creep on a vertical strike-slip fault
+\citep{Savage:Prescott:1978} and supershear dynamic spontaneous
+rupture of a 60 degree dipping normal fault in a Drucker-Prager
+elastoplastic medium. This second benchmark corresponds to benchmark
+TPV13 in the SCEC suite of dynamic spontaneous rupture benchmarks
+\citep{Harris:etal:SRL:2011}.
\subsection{Quasi-static Benchmark}
\label{sec:verification:quasi-static}
@@ -1408,7 +1413,9 @@ results against the analytical solution
results against the analytical solution of
\citet{Savage:Prescott:1978}. This problem consists of an infinitely
long strike-slip fault in an elastic layer overlying a Maxwell
-viscoelastic half-space. Figure~\ref{fig:savage:prescott::solution}
+viscoelastic half-space. The parameter files for this benchmark are
+available in the quasistatic/sceccrustdeform/savageprescott directory
+of the benchmark repository. Figure~\ref{fig:savage:prescott::solution}
illustrates the geometry of the problem with an exaggerated view of
the deformation during the tenth earthquake cycle. Between earthquakes
the upper portion of the fault is locked, while the lower portion
@@ -1519,7 +1526,9 @@ stress-drop, supershear, dip-slip earthq
(very large) ground motions, large slip, and fast slip rates
\citep{Harris:etal:SRL:2011}. It uses a Drucker-Prager elastoplastic
bulk rheology and a slip-weakening friction model in a depth-dependent
-initial stress field.
+initial stress field. The parameter files for this benchmark are
+available in the dynamic/scecdynrup/tpv210-2d and
+dynamic/scecdynrup/tpv210 directories of the benchmark repository.
Figure~\ref{fig:tpv13:geometry}
show the geometry of the benchmark and the size of the domain
@@ -2200,6 +2209,7 @@ VecMAXPY & 1 & 1.0 & 1819 \\
& Vp & 5.716 km/s \\
& Vs & 3.300 km/s \\
& Density ($\rho$) & 2700. kg/m$^3$ \\
+ & Nondimensional viscosity ($\eta^*$) & 0.4 \\
\hline
\end{tabular}
\tablenotetext{a}{Basic simulation parameters for the SCEC
diff -r 4c38d2c89757 -r ea734d8e36ef response_jgr.tex
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/response_jgr.tex Tue Feb 05 14:21:20 2013 -0800
@@ -0,0 +1,515 @@
+%-*- TeX -*-
+%
+% ----------------------------------------------------------------------
+%
+% Brad T. Aagaard
+% U.S. Geological Survey
+%
+% ----------------------------------------------------------------------
+%
+
+\documentclass{reviewresponse}
+
+% ==================================================================
+\begin{document}
+
+\maketitle
+
+% ----------------------------------------------------------------------
+\reviewer{Associate Editor}
+
+\comment{%
+ Reviewer \#1 requests clarifications on the stability criteria used,
+ also on the mathematical notation. These are important aspects that
+ need to be addressed.
+}{%
+ We added a discussion of the LBB stability criterion (details
+ discussed in Reviewer \#1 response). The AGU notation style uses
+ bold for tensors and matrices which is not as clear as our preferred
+ notation of using overbars. We adjusted the font markup in our LaTeX
+ file to ensure the Greek symbols are bold.
+}%
+
+\comment{%
+ Additionally, Reviewer \#2 remarks a short-coming regarding the use
+ of the Lagrange multipliers in the context of material discontinuity
+ across the fault (i.e. a bi-material interface). This aspect can be
+ accounted for / remedied, but the authors have to include this point
+ of concern/short-coming in their revised manuscript.
+}{%
+ The Lagrange multiplier implementation for fault slip does not
+ suffer the shortcoming claimed by Reviewer \#2. We added some
+ additional discussion (see details under Reviewer \#2 response) to
+ clarify how our implementation is similar to the often used
+ ``Traction at Split Node'' (TSN) implementation and differs from a
+ double couple point source implementation.
+}%
+
+\comment{%
+ In addition to the above comments, I would like to ask how material
+ complexities in 3D are handled in this implementation, and if
+ seismic wave propagation is included as well. Since the authors have
+ conducted a set of the 3D SCEC benchmark tests for dynamic ruptures,
+ I assume that 3D velocity-density models can be naturally included,
+ but it is not stated explicitly.
+}{%
+ Added a statement about how 3-D variations in physical properties are
+ handled at line ??.
+}%
+
+\comment{%
+ Additionally, seismic waves seem to be a natural by-product of the
+ dynamic calculations carried out, but it is of course a matter of
+ the domain size if these are computed out to distances of
+ seismological (earthquake-engineering) distances. A few comments
+ along these lines would be helpful
+}{%
+ Added a statement about suitability of the code for ground motion
+ simulations at line ??.
+}%
+
+\comment{%
+ L 40: not all models are "simplistic"; several researchers explored
+ random initial stress, variable Dc and other dynamic rupture
+ parameters (e.g. Ripperger et al, 2006; Ampuero et al, 2006; Dalguer
+ and Mai, 2011, among others)
+}{%
+ Removed ``simplistic'' from statement. Intended meaning is now clear.
+}%
+
+\comment{%
+ L 43: see Hillers et al, 2006, 2007 for earthquake-cycle simulations
+ with variability in rate-and-state parameters
+}{%
+ Added reference to Hillers et al. 2006 at line ??.
+}%
+
+\comment{%
+ L 106+112: which boundary conditions? not specified explicitly so far
+}{%
+ Explicitly mention which boundary conditions at line ??.
+}%
+
+\comment{%
+ L 113: define the double-dot product, " : "
+}{%
+ Clarified meaning of colon symbol.
+}%
+
+\comment{%
+ L 432: does the inverse of K exist for all cases?; the part of the
+ sentence with "K . L" remains unclear ...
+}{%
+ For this case of quasi-static simulations K is the stiffness matrix,
+ so the diagonal terms are all positive and their inverse exists. The
+ awkward sentence structure and notation was reworked to improve
+ readability.
+}%
+
+\comment{%
+ L 588: "the solve" --> replace with "the solution"
+}{%
+ Changed ``the solve'' to ``in solving the linear system of
+ equations``.
+}%
+
+\comment{%
+ Multiple times, the authors may remove "in order" in language
+ constructs like "in order to" ...for brevity
+}{%
+ Yes, the phrase ``in order to'' was used excessively.
+}%
+
+
+
+% ----------------------------------------------------------------------
+\reviewer{Reviewer \#1}
+
+\comment{%
+ The paper properly describes the implementation used. It has an
+ excellent description of the scalability and algorithmic details.
+ Many details on how to reproduce the simulation results are lacking,
+ but this can be easily fixed, particularly if pointers to the
+ particular problem's set ups are given to the PyLith repository.
+}{%
+ Added statement that the mesh generation and simulation parameter
+ files are available in the CIG subversion repository at line ??. The
+ directories for each the code verification benchmarks are also
+ added to the text.
+}%
+
+\comment{%
+ The major concern I have about the manuscript is that there is no
+ discussion of the stability conditions that the Lagrange multiplier
+ and bulk deformation discrete spaces need to satisfy. The authors
+ have omitted the discussion of the LBB conditions of the discrete
+ spaces (from the mathematical point of view) or a degrees of freedom
+ count (from the engineering point of view) to assess and ensure that
+ the discrete system is invertible. I believe the manuscript needs to
+ include a section where this space compatibility is discussed and how
+ is it controlled within PyLith. That is, in equations (130-(16) a
+ description of the rationale for the selection of $N_p$, $N_m$, and $N_n$
+ should be provided. For this reason is that I request a major
+ revision. Otherwise, the paper is well written and very descriptive.
+}{%
+ Added a discussion of the selection of basis functions for the
+ Lagrange multipliers and how they satisfy the LBB stability
+ criterion at lines ??--??.
+}%
+
+\comment{%
+ In equation (6) define the double contraction and dot product used.
+}{%
+ The double contractions is the double inner product.
+}%
+
+\comment{%
+ Paragraph below equation (6). Please define the tractions $l_{+}$ and $l_{-}$.
+}{%
+ Added description of the fault tractions.
+}%
+
+\comment{%
+ Does equation (8) include non-tangential slip? I believe your system
+ allows for crack openings. how do you account for them in your
+ formulation?
+}{%
+ Added statement that the constraint equation applies to both
+ tangential slip and fault opening.
+}%
+
+\comment{%
+ The notation used in equations (17) and (18) is very confusing to me.
+ That is, I know how the discrete system will look, but the notation
+ is not helpful. For example, in the first term of equation (17), what
+ is the transpose of the spatial gradient of the vector that contains
+ all basis functions? how do you take a double contraction of it with
+ the stress state? Similar concerns apply to each term of both
+ equations. In particular the third term of (17) and the first of
+ (18), what is the dot product of a vector with another vector, my
+ understanding is a scalar. Then what it means to take a dot product
+ of that scalar with a vector (transposed)? I recommend the authors to
+ use component (index) notation, such as the one used in Hughes's
+ finite element book, for example. Similar comments correspond to
+ equations (19), (20), (22) and (23).
+}{%
+ Improved font markup to ensure bold Greek symbols for tensors
+ (following AGU style guide). We prefer vector notation to index
+ notation and believe the current notation is clear.
+}%
+
+\comment{%
+ In page 9, just before section 2.1 The system's Jacobian (A) is
+ used without being defined. Additionally, in the last equation of
+ that paragraph, in general, A on the left-hand side does not need to
+ be the same one on the right-hand side of the equation, thus please
+ be more precise on notation usage.
+}{%
+ The purpose of this paragraph is to outline the general form for the
+ linear system of equations. In the following sections the details of
+ the terms are clearly described. We have adjusted the wording to
+ reflect the focus of the paragraph.
+}%
+
+\comment{%
+ Paragraph below equation (20), last sentence. What guesses are used
+ to compute the residuals in equations (19)-(20)?
+}{%
+ Added statement that the initial guesses for the solution are zero.
+}%
+
+\comment{%
+ Equation (21), please define what is the action of gradient
+ transposed. Additionally, this stress increment is not frame
+ invariant, thus constrains the theory to small deformations and
+ rotations.
+}{%
+ Defining strain in terms of the gradient and gradient transposed is
+ quite common and well-defined mathematically.
+}%
+
+\comment{%
+ Equation (34), what is the value of $\eta^*$ used and is $\frac{\partial
+ u}{\partial t}$ evaluated explicitly or implicitly? Please give the
+ value of $\eta$ for every simulation discussed in the paper.
+}{%
+ Added nondimensional viscosity parameter to Table 7. These are the
+ only simulations with explicit time stepping and the nondimensional
+ viscosity.
+}%
+
+
+\comment{%
+ Page 14. I have the following concern, how stable is the link to
+ Aagaard et al, 2012?? The paper under review will be published in
+ JGR, will the link to the PyLith user manual be available in the
+ future, or should the manual be published as supplementary material
+ by JGR??
+}{%
+ The link to the PyLith manual is stable. The PyLith webpage has had
+ the same URL for five years. All versions of PyLith released starting
+ with v1.0 in 2007 along with the corresponding manuals remain
+ available at the same URLs.
+}%
+
+\comment{%
+ Second paragraph in page 20, last sentence. The no slip at the fault
+ end is imposed weakly using Lagrange multipliers, thus do you observe
+ dependence of the simulation results on the ordering?
+}{%
+ Fixed wording to clarify what we mean by setting the slip to zero at
+ the buried edges of the fault.
+}%
+
+\comment{%
+ Equation (46) leads to an error of order O($\Delta t^2$), which for
+ explicit methods should be fine, since the lead to small time steps,
+ but further down the paper, large time steps are reported, how are
+ these two aspects reconciled?
+}{%
+ Added clarification that in quasi-static simulations the time
+ stepping is simply a series of static simulations at line ??.
+ }%
+
+\comment{%
+ Please try to more consistently use process, processor and core, it
+ took me a while to realize that these mean the same in the context
+ of this paper.
+}{%
+ Adjusted text to consistently refer to process, core, and processor
+ appropriately.
+}%
+
+\comment{%
+ Second paragraph in page 28, last sentence. What is the cost if you
+ include quadrature of residual, preconditioner formation and
+ application? Please discuss further the data shown in the top of
+ figure 6, I believe it is important.
+}{%
+ The cost of the numerical integration and preconditioner setup and
+ application are included in Figure 6. Added a statement clarifying
+ this to the caption of Figure 6.
+}%
+
+\comment{%
+ Page 20, last sentence. Why do you only consider iteration counts?
+ Set up and application can obscure this data. That is, if I apply
+ the LU factorization the iteration count is always 1, but it is too
+ expensive to form and apply.
+}{%
+ The discussion is focused on preconditioners we have found useful in
+ solving small problems. Their performance in terms of number of
+ iterations for large problems differs considerably, so we focus on
+ that metric before considering runtime in the parallel scalability
+ section.
+ }%
+
+\comment{%
+ Line 710, please explain how can the simulation be stable with a
+ time step size of 5 years, it seems to me that the explicit
+ stability limit without modification for Earth materials should in
+ the order of seconds to minutes, years is too many orders of
+ magnitude away? Please explain how are the material coefficients
+ modified to achieve this time step sizes.
+}{%
+ As described in the text, quasi-static simulations are a series of
+ static simulations. The accuracy of the solution is controlled by
+ the ability to resolve the temporal variation in the physical
+ properties and boundary conditions. In this case, the Maxwell time
+ controls the accuracy of the solution, so we discuss the time step
+ in terms of the Maxwell time.
+}%
+
+\comment{%
+ Page 37, lines 758-760. Effective pressures need to account for
+ drainage speed (pressure build up due to undrained conditions under
+ shaking). how are these affects accounted for?
+}{%
+ The SCEC Dynamic Spontaneous Rupture benchmark problem does not
+ account for drainage speed. The benchmark specifications are
+ extremely detailed, so we refer the reader to the benchmark
+ specification.
+}%
+
+
+% ----------------------------------------------------------------------
+\reviewer{Reviewer \#2: Sylvain Barbot}
+
+\comment{%
+ The first paragraph needs several representative citations. The work
+ of Ampuero (Ampuero \& Rubin "Earthquake nucleation on rate and
+ state faults - aging and slip laws" 2008), Matsuzawa (T. Matsuzawa
+ and H. Hirose and B. Shibazaki and K. Obara, "Modeling short- and
+ long-term slow slip events in the seismic cycles of large subduction
+ earthquakes" 2010), and Lapusta (Lapusta \& Barbot "Models of
+ earthquakes and aseismic slip based on laboratory-derived rate and
+ state friction laws" 2012), for example, should be acknowledged.
+}{%
+ ??
+}%
+
+\comment{%
+ The statement "influence earthquake rupture propagation, and the
+ dynamics of rupture propagation, in turn, affect postseismic
+ deformation." needs a citation. The paper by Chen \& Lapusta
+ ("Scaling of small repeating earthquakes explained by interaction of
+ seismic and aseismic slip in a rate and state fault model" 2009) is
+ a good illustration of the interaction between seismic and aseismic
+ processes. Observations of such processes can be found in the work
+ of Igarashi (T. Igarashi and T. Matsuzawa and A. Hasegawa "Repeating
+ earthquakes and interplate aseismic slip in the northeastern Japan
+ subduction zone", 2003), Ito (Y. Ito and K. Obara and K. Shiomi and
+ S. Sekine and H. Hirose "Slow earthquakes coincident with episodic
+ tremors and slow slip events", 2007) and many others really...
+}{%
+ ??
+}%
+
+\comment{%
+ A good illustration of models that reproduce multiple earthquake
+ cycles is the work of Barbot (S. Barbot, N. Lapusta and J.-P. Avouac
+ "Under the Hood of the Earthquake Machine: Towards Predictive
+ Modeling of the Seismic Cycle" Science 2012). This study involves
+ modeling of fault-induced deformation at the time scales going from
+ microseconds to several centuries and is the most advanced model
+ attempting to reconcile as many observations along the earthquake
+ cycle (geodesy along co-, post-, and inter-seismic periods,
+ seismicity, paleo-seismicity, precise hypocenter locations). The
+ model uses a persistent model of friction properties to explain
+ simultaneously many different aspects of the Parkfield earthquake
+ cycle. In the long run, we will be able to do similar studies with
+ Pylith and incorporate more realistic representations of the bulk
+ properties of rocks. This is an exciting prospect.
+}{%
+ ??
+}%
+
+
+\comment{%
+ The statement about the work of Kaneko et al [2011] (paragraph 3)
+ being some of the most advanced model seems partial. Could you
+ include a more inclusive list of advanced models of the earthquake
+ cycle. Of course Barbot et al. 2012 comes to mind, but you could also
+ include the work of Mai, or Dalguer, for example (A Gabriel, J.-P.
+ Ampuero, L. Dalguer, M. Mai "Source properties of dynamic rupture
+ pulses with off-fault plasticity" 2012) and from other people in the
+ broader community (not only US-centric!).
+}{%
+ ??
+}%
+
+\comment{%
+ A balanced citation of the community is important because we can
+ anticipate this paper to be quite influential, so it needs to lead
+ by example.
+}{%
+ We appreciate Sylvain's help in improving the choice of references
+ cited in the manuscript.
+}%
+
+\comment{%
+ In the formulation developed in Section 2, the Lagrange multipliers
+ for fault slip correspond to equivalent surface tractions. It is
+ clear from the double-couple theory (see Aki and Richards or Paul
+ Segall's book) that equivalent body force and equivalent tractions
+ scale with the elastic moduli. The displacements $u_i$ along a fault
+ with normal vector $n_j$ give rise to the potency tensor
+
+ \[
+ E_{ij} = (u_i n_j + u_j n_i) / 2
+ \]
+
+ and the equivalent body force
+
+ \[
+ m_{ij} = C_{ijkl} E_{kl}
+ \]
+
+ Once projected on the fault surface, this give the tractions
+
+ \[
+ t_i = m_{ij} n_j
+ \]
+
+ If the elastic moduli are different on each side of the fault, then
+ the equivalent tractions are different on each side. Assuming that
+ the Lagrange multipliers on each side of the faults, l+ and l-, are
+ in general equal and opposite is therefore incorrect.
+
+ This result is mathematically well established, but it can still be
+ confusing because it may appear that resulting tractions on each
+ sides of the fault could be different, which is against fundamental
+ physical insight from Newton's third law of dynamics stating that
+ for every action there is an equal and opposite reaction. But the
+ equivalent traction - the Lagrange multiplier - is not the final
+ traction on the fault resulting from solving the full dynamic
+ system. The equivalent traction is the force required to move the
+ fault by a prescribed amount. In the case of a bi-material
+ interface, prescribing different equivalent tractions on the
+ different sides of the fault is the way to obtain the same traction
+ at the end of the computation, and therefore satisfy Newton's law.
+
+ These considerations could benefit a benchmark calculation using
+ bi-material interface. The work of Barbot shows simple
+ two-dimensional benchmarks where a solver using equivalent tractions
+ is tested successfully against analytical solutions with dislocation
+ conditions (Barbot S., Fialko Y., Sandwell D., "Effect of a
+ compliant fault zone on the inferred earthquake slip distribution",
+ J. Geophys. Res., 113, B06404, 2008, doi:10.1029/2007JB005256).
+
+ The theoretical background presented in Section 2 should reflect
+ these considerations and the constraints of equal and opposite
+ Lagrange multipliers may have to be relaxed. This should not affect
+ most of the paper (bi-material contacts are not tested), but perhaps
+ performance.
+}{%
+ Sylvain's assertion that there is a flaw in the Lagrange multiplier
+ approach is incorrect. He is confused by the difference between
+ using the domain decomposition approach which uses the fault
+ traction as the Lagrange multipliers and the double couple point
+ source approach with an effective plastic strain. The tractions
+ resolved from body forces associated with the effective plastic
+ strain are not equal to the fault tractions; the tractions
+ coming from elasticity must also be included. These are included in
+ our formulation.
+
+ To clarify this distinction, we added a paragraph (lines ??--??)
+ discussing the similarities and differences between our domain
+ decomposition approach using Lagrange multipliers and other methods
+ of implementing fault slip, including ``Traction at Split Nodes''
+ and double couple point sources.
+}%
+
+\comment{%
+ A global search for "a a " or "the the " could fix some repetitions.
+}{%
+ Fixed.
+}%
+
+\comment{%
+ In Section 2.5, the description of the solver comes before
+ introducing what rheology is assumed. I think linear slip weakening
+ is assumed, but I'm not entirely sure.
+}{%
+ Added a statement that the solver is independent of the fault
+ constitutive model at line ??.
+}%
+
+\comment{%
+ In Section 2.5, you mention rheologies such as linear slip weakening
+ and Dieterich-Ruina rate-and-state friction. In the former, velocity
+ is unspecified and is a result of the optimization; on the later,
+ velocity is a simple function of stress and state variables.
+ Intuitively, it seems that different solvers could be better adjusted
+ to these different cases, but it's not discussed in this paragraph.
+}{%
+ We added a note that the solver could be
+ tuned to yield faster convergence for specific fault constitutive
+ models at line ??.
+}%
+
+
+% ==================================================================
+\end{document}
+
+% End of file
diff -r 4c38d2c89757 -r ea734d8e36ef response_usgs.tex
--- a/response_usgs.tex Wed Jan 23 07:41:39 2013 -0800
+++ b/response_usgs.tex Tue Feb 05 14:21:20 2013 -0800
@@ -4,8 +4,6 @@
%
% Brad T. Aagaard
% U.S. Geological Survey
-%
-% <LicenseText>
%
% ----------------------------------------------------------------------
%
More information about the CIG-COMMITS
mailing list