Computational Infrastructure for Geodynamics Wiki

1. Xianjie zha: where can I find the slides?

Brad Aagaard: From the workshop web page (https://geodynamics.org/cig/events/calendar/2015-cdm-tutorial). Click on “Agenda”. See the links labeled “slides” in the Title column in Session I.

18. Wajahat Ali: Thanks for your guidence I'm now running pylith on Linux its more conviniant on linux platform. Do you have any idea when could we access the recorded tutorials ?

Brad Aagaard: There are links to streaming versions of the recorded sessions on the agenda. Versions partitioned into smaller chunks and downloadable will be posted late next week (I have to play them back in real-time so it takes a while).

25. Justin R. Lindeman: Will next year's PyLith training be a workshop (e.g. 2014) or a tutorial (e.g. 2015)?

Brad Aagaard: We are planning an in-person workshop similar to previous versions (for example 2012, 2014) that combines training with science talks and discussions. We expect it to be the 3rd or 4th week of June 2016.

2. Stephanie Wollherr: what is the difference between the Lagrange multiplier approach and the traction-at-split-nodes approach for the slip implementation at the fault?

Brad Aagaard: In general there is very little difference in the mathematics between the Lagrange multiplier approach and traction-at-split-nodes. In the bookkeeping and implementation in actual code, there are differences. The cohesive cells and Lagrange multipliers provide a topological and mathematical description of the fault traction, which we find convenient.

3. what is the viscosity range that Pylith can compute ?

Brad Aagaard: PyLith includes nondimensionalization of all values internally, so selecting appropriate scales allows PyLith to efficiently solve problems across a vary wide range of scales. Domains from cm in size to thousands of km can be used. Similarly, one can use time scales from fractions of a second to thousands of years. As a result, PyLith can accommodate almost any reasonable viscosity.

20. Jeanne Sauber: How many processors/cores can they utilize?

Brad Aagaard: LU is serial, so it uses 1 core. MUMPS and SuperLU can use as many cores as you want. The main limitation is memory. Matt is suggesting to *start* with LU, but these are usually not the *best* solvers as Matt is discussing.

Memory bandwidth is generally the main bottleneck in finite-element computations, especially elliptic problems like quasi-static elasticity. As a result, using more compute nodes on a cluster is generally the best way to speed things up. Of course, if you spread things too thin then communication among processes becomes the bottleneck. Using more cores on a given machine often doesn't speed things up much because the memory bus is saturated.

24. Farrokh: Can we define the “minimum Jacobian” in PETSc? I ask this question, because I remember I got the error message that the Jacobian of my FE model was smaller than minimum Jacobian. It was a small numerical experiment for barzalian test (few inch of cylinder smaple), quasi static model, under velocity compressive boundary condition.

Brad Aagaard: The work Jacobian refers to two different things in PyLith. There is the Jacobian of the system of equations. This corresponds to the tangent stiffness matrix for an elasticity problem. There is also the Jacobian which maps a reference cell to the actual cell. The error you were getting is associated with the mapping. We check the determinant of the Jacobian against a threshold (default value is 1.0e-6). A value smaller than the threshold generally indicates something is wrong with the problem setup and would lead to very poor conditioning of the system of equations. Usually a small Jacobian means the length scale used in the nondimensionalization is too big (the default is 1.0 km).

8. so we can do forward modeling (if we have disp and x-y) ?

Brad Aagaard: I am not sure what you mean. The current version of PyLith is designed for forward calculations. It does have a user-friendly interface for static Green's functions that are intended for use in inverting for coseismic slip from GPS/InSAR data in complex 2-D and 3-D domains.

16. Is it possible to extract elastic medium parameters through the surface displacement using this code?

Brad Aagaard: PyLith does not currently have support for adjoint simulations for constraining elastic properties. This could be done as a series of forward problems. This may be possible in future versions of PyLith via the multiphysics capabilities that will make it possible to supply your own governing equation.

31. Farrokh: What do suggest to switch the analysis from implicit to explicit? If we need dynamic rupture analysis only in very small part of the whole time duration, and the most part it should be quasi-static, what do you suggest to handle it?

Charles Williams: Use quasi-static to get to point of spontaneous rupture, then use those as initial stresses for dynamic solution. We are working towards earthquake cycle modeling that would allow automatically switch back and forth between dynamic and quasi-static solutions. The first step towards that goal is better adaptive time stepping, which is scheduled for v3.0.

32. Wajahat Ali: Does any part of pylith is covering induce seismicity ?

Charles Williams: Upcoming multiphysics versions of PyLith will address this more completely (e.g., fully coupled poroelasticity). At present, you can estimate changes in effective normal stress and apply these as initial conditions to existing faults.

33. Wajahat Ali: if we don't have a known fault in an area can we use elastic properties from sonic logs to identify faults?

Brad Aagaard: PyLith does not contain any functionality to identify faults from sonic logs.

9. Wajahat Ali: Hi- I'm in pylith-2.1.0-linux-x86 folder can you please guide me how to proceed next to run Pylith through terminal ?

Brad Aagaard: cd src/pylith-2.1.0/examples. The step03.cfg that I demonstrated is in 3d/hex8. See the PyLith manual for where the folders associated with the other examples.

10. Wajahat Ali: so on linux system I first need to install using INSTALL file.

Brad Aagaard: We strongly recommend that you start with the PyLith binary. There are detailed instructions in the manual. The basic steps are: (1) unpack the binary. (2) Change to the top-level directory of the unpacked distribution, (3) run “source setup.sh”. You can then change to a directory with the example input files, such as examples/3d/hex8, and run a simulation using PyLith.

26. can we use echo— in script ?

Charles Williams: You cannot put shell or Python commands in the .cfg file. If necessary you could edit the Python code, but would then need to recompile it. We recommend using the Python and C/C++ debuggers rather than using print statements for debugging. The PyLith binary contains optimized code, so for debugging it is best to build from source.

4. what about easy restarting (i.e. if you run an EQ cycle model for n-years and want to continue for longer total time)?

Brad Aagaard: PyLith has good support for restarting for prescribed slip simulations via initial stresses and state variables. Spontaneous rupture (fault friction) has some ability to specify initial states. We don't have a mechanism for specifying initial displacements/velocities. We plan to rectify this in PyLith v3.0 as part of switching to PETSc time stepping.

We did find a bug in the initial stress specification for viscoelastic materials in v2.1 (I think it also applies to earlier releases too). This will be fixed in the upcoming release in the next few weeks.

12. Wajahat Ali: can you please explain more about nearest neighbor and will this step introduce error ? in locations

Brad Aagaard: For spatial databases the default is nearest point interpolation (interpolation is a misnomer because it isn't interpolation but just assigning values). The linear interpolation is conventional linear interpolation. Which one you uses depends on the type of spatial variation you have (linear, piecewise uniform, etc). You can optimize the choice and spatial database for the specific variation that you want.

14. Farrokh: Can we use the same procedure you described for multiple fractures? For example is we have multiple opening fractures under the ground, and by using surface deformation to find opening distributions along multiple parallel fractures

Brad Aagaard: Yes, you can use the Green's functions procedures for inversions on multiple faults. In a single PyLith simulation you can only compute Green's functions for a single fault, so multiple simulations would be required. As a result, it involves more bookkeeping than a single fault, but the procedure is generally the same.

15. Katie Jacobs: Are there any general rules about how many green's functions we would want to generate over a given distance?

Brad Aagaard: This really depends on the number of observations and the spacing between observations. Checkerboard resolution tests can give you a general idea.

28. Farrokh: My question is about modeling the dynamic rupture considering insitu stress. Since for the explicit dynamic analysis Pylith use a trivial solver, I was wondering how can include static calculation before time-stepping? Can we use db_initial_stress to include initial stress throughout the whole domain? I was using db_initial_stress applied along the boundaries to include insitu remote stresses; however, the results I got showed applying loads on the boundaries generate a heavy P-wave that made the results wrong.

Brad Aagaard: For a linear elastic model the deformation depends only on the change in stresses not the initial stresses so we usually don't use them. For simulations with spontaneous rupture (FaultCohesiveDyn), we need the initial tractions on the fault, so we usually impose them directly if we have a linearly elastic material. For elastoplastic and viscoelastic models the initial stresses in the bulk medium can be important. For dynamic spontaneous rupture simulations, we don't tend to use viscoelastic materials; we do use initial stresses for elastoplastic problems with the Drucker Prager bulk consitutive model.

29. Farrokh: There were some additional settings specific to fault friction. in [pylithapp.petsc] for dike intrusion (step20.cfg). I was wondering if we want to include an opening mode pressurized fracture in a dynamic rupture (explicit time stepping), since we don’t use PETSc, do we still need to use these modifications for fault friction?

Brad Aagaard: The zero_tolerance parameter is relevant for explicit time stepping, but the PETSc solver settings are not because we don't use the PETSc solvers with explicit time stepping. For dike intrusions the open_free_surface parameter would also be relevant.

30. Farrokh: Can we use nonuinform, automatic time step in spontaneous dynamic rupture?

Brad Aagaard: PyLith uses the central difference scheme for explicit time stepping. The stable time step is controlled by the CFL condition (time it takes the P wave to propagate across a cell), so a uniform time step is used. When we switch to using PETSc time stepping algorithms (planned for v3.0) we will have more options for adaptive time stepping in quasi-static simulations.

11. Wajahat Ali: sorry just to clear my concept slip along fault and displacements (x-y-z) are different?

Brad Aagaard: In 2-D the HDF5 output only includes the two components of motion (x and y components for displacement/velocity and left-lateral and opening for slip/tractions). ParaView automatically calls these two components x and y. Within ParaView we often want a vector (3 components) so we use the Calculator tool to assemble a vector from the components (displacement_x*iHat + displacement_y*jHat).

13. Farhan Javed 2: is it possible in pylith to output the displacement vector along the line of sight (LOS)?

Brad Aagaard: PyLith outputs displacement vectors (3 components in 3-D), so to get the LOS vector you can dot the displacement vector in output with the look direction.

5. How did you realize which type of mesh to choose ?

Brad Aagaard: We often decide on the type of mesh (quad vs tri and hex vs tet) based on the geometry. Hex/quad meshes are good for regular sized domains with relatively uniform discretization sizes. Tet/tri meshes are good for irregular domains and large variations in discretization sizes. In general you need to know something about the geometry of your domain and the spatial variation of the expected solution.

6. Could you explain more about how to deal with topography by UV ?

Brad Aagaard: We don't have time to cover complex meshing in this tutorial. There are some included with PyLith in examples/meshing. For topography see examples/meshing/surface_nurbs/dem. These examples are discussed in previous tutorials with videos on demand. See the CDM2014 tutorial (Mon 1:00pm-2:30pm) CUBIT/Trelis: Complex geometry and sizing functions at https://geodynamics.org/cig/events/calendar/2014-cdm-workshop/meeting-info/agenda/ and CDM2013 Session II: Intermediate CUBIT Meshing Strategies at https://wiki.geodynamics.org/software:pylith:cdm2013.

7. In Trelis, to model multiple fractures in a 3D body, do we need to define multiple volume including or Trelis can understand fractures as a surface of discontinuity in a single volume without putting fractures between adjacent volumes?

Brad Aagaard: PyLith relies on an interior surface (group of vertices associated with cell faces) for faults. CUBIT/Trelis can only generate this type of information if the surfaces for the faults (fracture surfaces) separate two volumes. This means the surfaces must extend to the edges of the volumes. Sometimes it is convenient to insert artificial surfaces (for example, a horizontal surface near the top of the mantle) so that you can divide your domain into pieces and allow your fault surfaces to truncate at this artificial surface rather than extend them all the way to the bottom of your domain. In either case you need not define the nodeset in CUBIT/Trelis to be the entire surface, but can create a set of vertices over a subset of the surface to match the extend of the fault/fracture.

19. Renier Ladron de Guevara: I've been trying to reproduce the volcano with dike and magma chamber with a 90m SRTM DEM model amd I am getting the errors I am copuing below. I will appreciate any help regarding the possible reason for this issue. The error comes out when running the dem2lines.py script.

Traceback (most recent call last):  File "./dem2lines.py", line 473, in <module>
    app.run()  File "/Users/Renier/Documents/PhD/Pylith_tutorial/pylith-2.1.0-darwin-10.6.8/lib/python2.7/site-packages/pythia-0.8.1.15-py2.7.egg/pyre/applications/Application.py", line 42, in run
    shell.run(*args, **kwds)  File "/Users/Renier/Documents/PhD/Pylith_tutorial/pylith-2.1.0-darwin-10.6.8/lib/python2.7/site-packages/pythia-0.8.1.15-py2.7.egg/pyre/applications/Shell.py", line 143, in run
    method(*args, **kwds)  File "/Users/Renier/Documents/PhD/Pylith_tutorial/pylith-2.1.0-darwin-10.6.8/lib/python2.7/site-packages/pythia-0.8.1.15-py2.7.egg/pyre/applications/Stager.py", line 19, in execute
    return self.main(*args, **kwds)  File "./dem2lines.py", line 120, in main
    self._resampleDem()  File "./dem2lines.py", line 266, in _resampleDem
    zTmp = numpy.take(self.zIn, yIndices, axis=0)  File "/Users/Renier/Documents/PhD/Pylith_tutorial/pylith-2.1.0-darwin-10.6.8/lib/python2.7/site-packages/numpy/core/fromnumeric.py",
    File "/Users/Renier/Documents/PhD/Pylith_tutorial/pylith-2.1.0-darwin-10.6.8/lib/python2.7/site-packages/numpy/core/fromnumeric.py", line 103, in take
    return _wrapit(a, 'take', indices, axis, out, mode)  File "/Users/Renier/Documents/PhD/Pylith_tutorial/pylith-2.1.0-darwin-10.6.8/lib/python2.7/site-packages/numpy/core/fromnumeric.py", line 38, in _wrapit
    result = getattr(asarray(obj),method)(*args, **kwds)
    IndexError: index 2 is out of bounds for axis 0 with size 1

Brad Aagaard: This looks like a problem with getting the data into the arrays. The array sizes are different than expected.