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).

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.

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.

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.

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.

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.

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.