Hello,
I am trying to run a simple post-seismic relaxation model that I created with power-law viscoelastic lower crust/mantle. Below are the first few lines in my .cfg file (step02_postseis_powerlaw.cfg) where I define the time step:
[pylithapp]
dump_parameters.filename = output/step02-postseis-powerlaw.json
problem.progress_monitor.filename = output/step02-postseis-powerlaw-progress.txt
----------------------------------------------------------------------
problem
----------------------------------------------------------------------
[pylithapp.problem]
For this problem we must switch to a nonlinear solver.
implicit.solver = pylith.problems.SolverNonlinear
[pylithapp.problem.formulation]
time_step = pylith.problems.TimeStepAdapt ; Change the time step algorithm
[pylithapp.problem.formulation.time_step]
total_time = 20.0year
max_dt = 1.0year
adapt_skip = 10 ; Default value
stability_factor = 2.0 ; Default value
However, when executing the scripts (pylith step02_mat_powerlaw.cfg step02_postseis_powerlaw.cfg solver_fieldsplit.cfg), it appears that the time step is way too small compared to the 20-year total time (in the order of second):
/Users/Alpine/PyLith/pylith-2.2.2-darwin-10.13.6/lib/python2.7/site-packages/pylith/problems/TimeDependent.py:195:run
– timedependent(info)
– Preparing to advance solution from time t=0s to t=2.75653s.
I wonder why pylith is assigning such a small time step and how can I increase it (to maybe 0.5 year)?
I have tried using a uniform time step but it ended up in errors:
RuntimeError: Current nondimensionalized time step of 5.0000e-03 exceeds the nondimensionalized stable time step of 1.7470e-09.
application called MPI_Abort(MPI_COMM_WORLD, -1) - process 0
Thanks in advance for any help,
-Guo Cheng
The time step size is dependent on the relaxation time for the material. The relaxation time for a power-law material is described in the manual, and is dependent on both the material properties and the stress state. I suspect that your effective viscosity is very low, which will limit the size of time step that you can take. It also appears that maybe you should be using a smaller value for nondimensionalizing time.
Depending on the problem, the effective viscosity will generally increase with time, as the stresses relax, which will allow longer time step sizes as the simulation runs.
Cheers,
Charles
1 Like
Hi Charles,
Thank you for your reply. By “using a smaller value for nondimensionalizing time”, do you mean a smaller number for total_time, which should not be too large compared to the relaxation time? For example, I have also tried a simple Maxwell model with viscosity = 1e17, which has a relaxation time of ~7 days. In this case I probably should not set the total time to more than a year, correct?
If I do want to test out the long-term response, say 10 years, is it possible to run Pylith in parallel to increase the computing speed?
Best,
Guo
Nondimensionalization is discussed in section 4.1.3.1 of the Pylith manual (v2.2.2). The time scale used for nondimensionaliation in quasistatic problems should be on the order of the relaxation time.
You can run PyLith in parallel. On a multicore laptop or desktop machine you can often get about 2-3x speedup. How much speedup you get depends on the problem size, memory architecture, and number of cores.
It also appears to me that your relaxation time is much too small for a realistic geological problem. What sort of power-law material properties are you using?
Attached figure shows my test model, which consists of 3 domains: an elastic upper crust (blue), a viscoelastic lower crust/mantle (yellow), and an elastic slab (green). I have a fault plane with prescribed slip distribution in within the upper crust (highlighted in orange within the magenta domain).
For the power-law properties of the lower crust/mantle, I simply modified the power-law properties given by the example 7.9.6.6. Below is the .spatialdb file generated after running powerlaw_gendb.py. For the elastic properties, I use density, Vs, Vp of 3300 kg/m3, 4403.9 m/s, and 7627.7 m/s.
#SPATIAL.ascii 1
SimpleDB {
num-values = 3
value-names = reference-stress reference-strain-rate power-law-exponent
value-units = Pa 1/s none
num-locs = 10
data-dim = 1
space-dim = 3
cs-data = cartesian {
to-meters = 1
space-dim = 3
}
}
0.000000e+00 0.000000e+00 0.000000e+00 1.609650e+19 1.000000e-06 1.500000e+00
0.000000e+00 0.000000e+00 -5.000000e+04 3.781188e+12 1.000000e-06 3.500000e+00
0.000000e+00 0.000000e+00 -1.000000e+05 5.798732e+08 1.000000e-06 3.500000e+00
0.000000e+00 0.000000e+00 -1.500000e+05 1.159815e+07 1.000000e-06 3.500000e+00
0.000000e+00 0.000000e+00 -1.900000e+05 1.821615e+06 1.000000e-06 3.500000e+00
0.000000e+00 0.000000e+00 -2.100000e+05 9.095997e+05 1.000000e-06 3.500000e+00
0.000000e+00 0.000000e+00 -2.500000e+05 3.053179e+05 1.000000e-06 3.500000e+00
0.000000e+00 0.000000e+00 -3.000000e+05 1.131358e+05 1.000000e-06 3.500000e+00
0.000000e+00 0.000000e+00 -3.500000e+05 5.442705e+04 1.000000e-06 3.500000e+00
0.000000e+00 0.000000e+00 -4.000000e+05 3.103695e+04 1.000000e-06 3.500000e+00
Hi Brad,
Is this relaxation time a value calculated from the material properties I assigned (e.g. viscosity/elastic modulus for Maxwell) or a rough number that I estimate for my earthquake case (e.g. I assume my coseismic strain will be fully relaxed in 100 years)?
Thanks,
Guo
The nondimensionalization time scale given as a relaxation time is a single scalar value you provide in the parameter file separately from material property values. The normalizer scales are used to nondimensionalize all dimensional input quantities.
The relaxation time for a power-law material is given by equation 5.82 in the manual. You can estimate it for the different parts of your model by looking at the elastic stress values. You would then need to compute J_2’ from the stress and put that into the equation. Are you using gravity in your problem? If so, just the differential stresses from gravity could give you a huge value for J_2’, which would give you a very small relaxation time.
I don’t think I have incorporate gravity in my problem. How would you calculate differential stress sigma_d in calculating J_2’ in a case without considering gravity? Is it the difference between the overburden and the confining pressure?
If you’re not using gravity you don’t have to worry about that. In that case, your stresses are due purely to the applied dislocation. You can just look at the elastic solution and figure out where the stresses are really high. If you have very high stresses in regions where you have a low reference stress, this will lead to extremely short relaxation times. If you think your viscoelastic properties are correct, you will just need to let the stresses relax, and the relaxation times will gradually increase. Otherwise, you might want to consider adjusting your properties.
1 Like
