Computational Analysis/Plan for GEM

Find this at http://www.npac.syr.edu/users/gcf/gemcomp.html

Statement of Problem (from a computational point of view)

• Let i label the basic degrees of freedom which correspond to fault segments with the independent variables corresponding to the slip vector s(i or x,t) or more precisely the slip deficit fi(i or x,t).
• The slip deficit represents the deviation from the average (global) movement of the fault.
• In the original Rundle 88 calculation, there N=80 fault segments which averaged 30 km(horizontal) by 18 km(vertical) in size. This corresponds to a total fault length of 2400km.
• One would expect to need a time and position varying segment size with 100 metres as a goal for size in critical regions. This would increase number of segments by a factor 60,000 to a total of 5 million
• The paper Dynamics of a Traveling Density Wave Model for Earthquakes describes problem in the most general form we would need initially. This is formulated as a set of Green's function equations for the stress and displacements at the fault segments.
• We can "think" about fault segments as "particles" interacting with a long range force. This force is a "double-couple" i.e. a pair of equal and opposite dipoles.
• This force is just the Green's function linking the positions corresponding to the fault segments.
• This "force" G/T(i,j) for stress and displacement is in principle a function of time but this effect can be ignored in our initial computations.
  • Thus an unusual feature (compared to molecular dyamics or astrophysical particle simulations) is that forces on segment i are linear in displacements fi(j,t) at site j and coefficient in linear formula is fixed in time.
  • However later calculations would need to include the radiative time dependence in G/T(i,j) at some level.
• The most general approach is to evolve a set of N ordinary differential equations -- one for each "particle" i -- where each right hand side has a sum of N terms corresponding to long range force structure.
  • At each time step, one must examine the stresses and compare with the friction law to see if a slip will occur.
  • There are many different proposals for the friction law and the computational system needs to be set up to be able to flexibly change and compare results from different friction laws.
• An "earthquake" is a correlated ( i-m .. i-m+1 ... i i+1 .. i+n-1 i+n all "slip") set of segments that slip together.
• Earthquakes have various sizes with a "large" earthquake having depth 20 km and length 400 km ( this could be a million correlated points)

Choices in Solving GEM Equations

• We can formulate this 100 --> 5 million particle problem as either a particle dynamics or Monte Carlo method.
  • Given the analogies to statistical physics problems explored by Klein and Rundle, the possible interest in Monte Carlo methods appears natural although it has not been explored upto now.
• There are two important simplifications that have been used in current work
  • In the "Slider-Block" approach one only keeps the "nearest-neighbor" interactions and this approximation can be generalized to a screened approximation where Green's function elements G/T(i,j) are neglected when i and j are separated by an appropriate distance. With such a screened force, the computation becomes O(N) of course and not the O(N^2) of traditional long range force case
  • In the "cellular-automata" approach (used in Rundle 88), one uses a rather large time step ( 5 years in Rundle 88) and evolves the "lattice" of fault segments in discrete steps.
• In the "cellular-automata" approach, one sees a formulation very similar to that of statistical mechanics with an earthquake in GEM corresponding to a long range correlation forming a cluster in Statistical physics case.
• Clustering problems are rather hard to parallelize (we have experience with this from Swendsen Wang and other statistical mechanics problems -- see for instance my book "Parallel Computing Works!")
  • There are often not enough points in the cluster (earthquake) to allow much parallelism. In GEM this depends on total number of clustered points at each time step -- it must be zero or a lot to allow efficient parallelization
  • As earthquakes (clusters) are not uniform in physical domain (at given time step), one must use some sort of cyclic (block scattered) decomposition so that each processor participates roughly equally in any given earthquake scenario.
  • The parallelism is hardest in case of a screened or O(N) force computation as in long range force version, one can get parallelism in force computation and more generally the computation of forces (which is easy to parallelise) becomes a larger part of the computation time.
• Consider the most general travelling wave formulation where we integrate the differential equations
  • In the corresponding astrophysical problems, the computations have been revolutionized by the so called fast multipole methods which reformulate problem as either a particle - multipole (O(NlogN) ) or multipole - multipole ( O(N) ) method
  • These new methods have allowed tens of millions of particles to be simulated which appears to a size of interest for GEM. Note that molecular dynamics for some Chemistry problems still use O(N^2) approach as fast multipole methods only start to outperform the brute force approach in the regime of 10,000 particles.
  • The multipole methods naturally implement a variable spatial resolution for GEM where we only need fine resolution where the action is.
  • The manifestation of "clustering" problem in the cellular automata approach, becomes in this general formuation a variable time step with smaller time steps needed in earthquake region.
  • There are clearly several non trivial load balancing issues which are produced by variable time and space resolutions and actions (is there an earthquake or not)

Proposed Actions in GEM Computational Effort

1. Develop a flexible object-oriented friction package
2. Develop an approach to conveniently specify geometery of faults and choice of variable resolution segments
3. Define visualization and data input requirements for GEM and choose appropriate implementation for these features which can be used in any proposed computational solution system
4. Integrate 1) and 2) into a conventional cellular automata package allowing flexible choice of force screening
5. Analyse computational complexity of GEM to define load balancing issues for both cellular automata and travelling wave approaches.
6. Develop parallel version of cellular automata solution system
7. Integrate 1) and 2) into a simple brute force Travelling Wave solution system with fixed (in time) spatial resolution (which can however vary with position in fault system) and fixed position independent time steps
8. Extend 6 for variable time and space resolutions
9. Develop a parallel version of brute force travelling wave problem which is easy except for variable resolution issues.
10. Analyse and Develop fast multipole algorithms for GEM exploiting time independence of G/T(i,j) and special "double-couple" structure.
11. Implement sequential fast multipole solver for GEM allowing variable time steps.
12. Implement parallel fast multipole solver for GEM allowing variable time steps.
13. Modify above for any experimental or theoretical discoveries which could perhaps imply need for time dependent Greens functions or other complexities or changes. (The disclaimer clause allowing all plans to be changed!)