[CIG-SHORT] long-term simulation of subduction seismic cycles with rate and state friction law

[Brad Aagaard]

Shaoyang,

Make your own time step object (for example, you might call it 
TimeStepAdaptSlipRate) and put it in a separate Python file. A good 
starting point is TimeStepAdapt.py in src/pylith/problems. Note that you 
*must* have the time_step function return your object in order for 
Pylith to use it. This would be the function that gets called when you 
use the following in your .cfg file:

time_step = pylith.problems.TimeStepAdaptSlipRate

In the above, I am assuming you would put the file in 
$PREFIX/lib/python$PYTHON_VERSION/site-packages/pylith/problems (PREFIX 
is where PyLith is installed and PYTHON_VERSION is the version number of 
Python used). A better place to put the Python file with your new object 
is in a separate directory outside the PyLith distribution, so that you 
can easily maintain it separately from PyLith. Add this directory to 
your PYTHONPATH environment variable and add a __init__.py file (it can 
be empty or see the ones in the PyLith distribution) to this directory.

In setting up an adaptive time step, you can access the portions of the 
output files that have already been written just like you can from a 
Python script after the simulation is over. You may want to limit your 
I/O using slices. The index -1 is used to get the last slice, so 
slip_rate[-1,:,:] gets the last time slice for the slip rate array. Here 
is an example:

h5 = h5py.File("output/mysim-fault.h5", 'r', driver='sec2')
slipRate = h5['vertex_fields/slip_rate'][-1,:,:]
h5.close()

Regards,
Brad


On 11/25/2013 10:21 AM, Shaoyang Li wrote:
> Dear Brad,
>
> thanks for your reply. I am going to implement variable time step in
> quasi-static regime as you suggested and follow part of a paper of Kaneko et al
> JGR 2011.
>
> I am not sure how to alter Pylith code now. But as I guess, it may work just
> write a new time step method based on the TimeStep.py and TimeStepUser.py. I
> find out there are two folds in Pylith containing these two python scripts
> (not complied in /src/pylith/pylith/problems and compiled in
> lib/python2.7/site-packages/pylith/problems). I am not sure in right track and
> wandering in which fold to put new time step method?
>
> Another question is how to refer slip rate on fault during simulation? I am only
> knowing to look it in the output files after the simulation is done. Thanks a lot.
>
> Best regards
> Shaoyang
>
>
>
> On Thu, 21 Nov 2013 17:55:29 -0800
>    Brad Aagaard <baagaard at usgs.gov> wrote:
>   > On 11/21/2013 08:04 AM, Shaoyang Li wrote:
>   >> Dear Brad,
>   >>
>   >> I am trying to manually implement "adaptive" time step
>   >>of rate and state
>   >> friction law with some Python scripts.
>   >>
>   >> Now in Pylith, if using implicit method and constant
>   >>time step (e.g. dt = 2
>   >> year), the earthquake can happen in several time step
>   >>(e.g. rupture in 4 steps).
>   >>
>   >> So my plan is: simulate several seismic cycles with
>   >>implicit method and constant
>   >> large time step. Choose one stable seismic cycle and
>   >>export all the state
>   >> variables of the model just several time steps before
>   >>the cycle starts. Import
>   >> all the variables and decrease the time step to restart
>   >>a new simulation.
>   >>
>   >> Now, my problem is: when I only import the stresses (of
>   >>the time step just
>   >> before an mature seismic cycle) to restart new
>   >>simulation, the model again need
>   >> a lot time to evolve to the first seismic cycle while
>   >>NOT continue as the
>   >> previous calculated cycles. So I am wandering the state
>   >>information (e.g.
>   >> velocity, state variable etc ) of RS friction fault
>   >>should also be import. But
>   >> 1, the state variable is not support in output fault
>   >>information.
>   >> 2, can I import velocity for every nodes on fault as
>   >>initial state?
>   >
>   > You can output the state variable using the fault output
>   >and set it using a spatial database for the initial
>   >state. You need to do this for any state variable in the
>   >bulk or fault constitutive models.
>   >
>   > Currently, there is no mechanism to specify the initial
>   >displacement or velocity for the entire domain or the
>   >slip or slip rate on the fault. This could be a problem
>   >unless you can get the solution to rapidly spin up in a
>   >few time steps.
>   >
>   > If you are using quasi-static simulations, the
>   >alternative is to use a variable time step. The time
>   >stepping object is all in python so you can add a new one
>   >if you have a simple algorithm (for example, looking at
>   >the slip rates in the fault output).
>   >
>   > Regards,
>   > Brad
>   >
>