Solving Nonlinear Partial Differential Equations Arising from Earthquake Models on Distributed Memory Parallel Machines
Computer Science 240A Applied Parallel Computing, Spring 2009
University of California Santa Barbara

Problem Description

In the late 1960s, Burridge and Knopoff (BK) introduced a model that exhibited some characteristics similar to the dynamics of an earthquake fault [2]. They were interested in the role that friction plays in regard to the earthquake mechanism and how successive events relate to each other in time and space. The basic configuration of the BK model consists of several blocks coupled by springs to a moving loader plate representing the other side of the fault (see Figure below).

The equations of motion for this model include a friction term accounting for the roughness of the surface upon which the block slips. In the early 80s, rate and state dependent frictions laws were shown to be qualified in reproducing some behavior similar to that of earthquake faulting. Burridge and Knopoff incorporated a simple linear friction term in their model that was dependent on the block's velocity. The crucial improvements to the BK friction law were made by Dieterich, Ruina, Rice and others based on empirical studies of rock friction in the laboratory. They discovered that a nonlinear friction term with the incorporation of a state variable enabled the model to exhibit almost entirely the observed seismic behaviors such as stick-slip phenomena, fault healing and memory effects (see [4] for more on this).

We study a one-dimensional chain of spring-connected blocks elastically coupled and driven by a plate moving at a constant rate V_0 . The blocks slide along a rough surface according to the nonlinear Dieterich-Ruina (DR) friction law. We take the continuum limit and obtain the following elastic wave equation involving the slip, u(t,x) , the slip velocity, u_t (t,x), and the state variable \theta :

\label{eqn1} \left. \begin{array}{l} u_{tt} = c^2 u_{xx} - \gamma ^2 \left( {u + \frac{1}{\xi }\left( {\theta + \log (u_t + 1)} \right)} \right) \\ \theta _t = - (u_t + 1)\left( {\theta + (1 + \varepsilon )\log (u_t + 1)} \right) \\ \end{array} \right\}

where c^2 is the square of the wave speed, \varepsilon = (B - A)/A ( A and B are frictional parameters in D-R friction) measures the sensitivity of the velocity relaxation, \xi is the nondimensional spring constant, and \gamma is the nondimensional frequency. See [3] for a better description of these parameters.

Numerical Analysis

Dieterich-Ruina style friction in the spring-block model involves a logarithmic term whose nonlinearity has introduced additional difficulty in solving the problem. Due to the nonlinear term, analytic integration cannot be done even in the simplest case when one block is used. And while numerical solutions can be computed, the logarithmic term still proves to be a great challenge.

As discussed in [3], the system (1) has large negative eigenvalues and so it is considered stiff. Typically, the necessitates the use of implicit methods in order to maintain stability due to sharp transients in the problem itself, which in turn require extremely small time steps in order to resolve the solution. The uses of such methods require additional work often by way of Newton Iterates that involves evaluations of Jacobian matrices which are generally expensive to compute and difficult to parallelize.

However, it was determined in [3] that the nonlinear logarithmic term forces small time steps as the velocity approaches -1, where the logarithmic terms goes to - \infty and thus an implicit method would not give us the benefit of taking larger time steps. So, we can use an explicit method to solve the system, which requires less work than implicit methods, however have smaller absolute regions of stability. Since we are forced to take small time steps already, we can choose an explicit method whose region of absolute stability contains the largest value step size allowed by the nonlinear term.

In the past, the Dieterich-Ruina friction term has been altered because of the difficulties imposed by the nonlinear term. Either this alteration of the nonlinear term or the lack of better solution algorithms may explain why chaotic regimes have rarely been observed in simulations done with this friction law and realistic parameter values. We plan to use the full nonlinear logarithmic term taking as small a time step as necessary to maintain accuracy. Our assumption is that the parallelization of this problem will reduce the time required for numerical integration.

If we discretize the spatial derivative in (1) according to the Method of Lines the original PDE reduces to a set of ordinary differential equations (ODE). In addition, the second order time derivative can be reduced to first order through substitution. The resulting system of first order ODEs (2) can then be solved using any class of numerical schemes which are well understood for these types of problems.

\label{eqn2} \left. \begin{array}{l} \dot u_i = v_i \\ \dot v_i = \frac{{c^2 }}{{\Delta x}}(u_{i + 1} - 2u_i + u_{i - 1} ) - \gamma ^2 \left( {u_i + \frac{1}{\xi }\left( {\theta _i + \log (v_i + 1)} \right)} \right) \\ \dot \theta _i = - (v_i + 1)\left( {\theta _i + (1 + \varepsilon )\log (v_i + 1)} \right) \\ \end{array} \right\}

In general, current ODE solvers are implemented to solve problems of the form:

\dot y(t) = f\left( {t,y(t)} \right),y(t_0 ) = y_0 ,y:\Re \to \Re ^n ,f:\Re \times \Re ^n \to \Re ^n

where the function f is often nonlinear. Regardless of the numerical scheme, the function will need to be evaluated at least once each iteration.

Two methods we are considering are explicit Runge Kutta (RK) Methods and Adams-Bashforth (AB). RK methods utilize multiple stages to increase the order of the method. However, for some desired order of accuracy there is as many function evaluations, up to order 4. Also, each of the stages may depend on one another. Alternatively, the AB Methods only require a single function evaluation per iteration, after some initial startup.

However, AB has smaller stability regions than the RK methods. In addition, AB Methods require additional storage for past function evaluations from which it attains higher order. This may be an issue in terms of collecting the necessary information for post processing analysis; we may run out of memory. See [1] for more details on these methods.

Sources of Parallelism

If we examine the expanded system for the velocity in the set of equations (2), we see that other than the matrix-vector multiplication, each of the spatial trajectories is independent to one another.

\frac{d}{{dt}}\left[ {\begin{array}{*{20}c} {v_1 } \\ {v_2 } \\ \cdot \\ \cdot \\ \cdot \\ {v_m } \\ \end{array}} \right] = \left[ {\begin{array}{*{20}c} \beta & \alpha & {} & {} & {} & {} \\ \alpha & \beta & \alpha & {} & {} & {} \\ {} & \cdot & \cdot & \cdot & {} & {} \\ {} & {} & \cdot & \cdot & \cdot & {} \\ {} & {} & {} & \cdot & \cdot & \cdot \\ {} & {} & {} & {} & \alpha & \beta \\ \end{array}} \right]\left[ {\begin{array}{*{20}c} {u_1 } \\ {u_2 } \\ \cdot \\ \cdot \\ \cdot \\ {u_m } \\ \end{array}} \right] + \left( {\frac{{ - \gamma ^2 }}{\xi }} \right)\left[ {\begin{array}{*{20}c} {\theta _1 + \log (v_1 + 1)} \\ {\theta _2 + \log (v_2 + 1)} \\ \cdot \\ \cdot \\ \cdot \\ {\theta _m + \log (v_m + 1)} \\ \end{array}} \right],\alpha = \frac{{c^2 }}{{\Delta x}},\beta = - \left( {\frac{{c^2 }}{{\Delta x}} + \gamma ^2 } \right)

Thus, we can distribute the spatial points across p processors and implement a parallel matrix-matrix vector multiplication method which takes advantage of the banded structure of the matrix. It is clear that the two remaining equations in (2), the state and substitution variables, also have independent spatial trajectories. So, we expect near p speed up in the function evaluation.

[6] and [5] discuss ways to parallelize the stages in RK. Thus, if we can distribute each stage on as many processors, we can reduce the time to compute the function evaluations to that of just one. In addition, we obtain the added benefit of better stability and lesser memory requirements. [7] discusses parallel issues in AB.

Objective

We will examine RK, AB, and other methods, to find suitable schemes with which we can parallelize well. Our goal is to parallelize (2) and obtain strong scaling up towards 32 processors. We will implement and test using Message Passing Interface (MPI) in C on either San Diego Supercomputing’s Teragrid or OnDemand. In addition, we will design with the ability to extend to higher spatial dimensions as needed.

References

[1] Ascher, U. M., Petzold, L. R.: Computer methods for ordinary differential equations and differential-algebraic equations, SIAM, Philadelphia, 1st edition, 1998.

[2] Burridge, R., and Knopoff, L: Model and theoretical seismicity, Bull. Seism. Soc. Am., 57, 341--371, 1967.

[3] Erickson, B., Birnir, B.:The Scalings of a Burridge-Knopoff Wave Equation with Detererich-Ruina Friction. ts. University of California Santa Barbara. 2009.

[4] Erickson, B., Birnir, B. and Lavallee, D.: A model for aperiodicity in earthquakes, Nonlin. Processes Geophys., 15, 1--12, 2008.

[5] Korch, M., Rauber, T.: Optimizing locality and scalability of embedded Runge-Kutta solvers using block-based pipelining. J. Parallel Distrib. Comput. 2006: 66(3): 444-468.

[6] Korch, M., Rauber, T.: Scalable Parallel RK Solvers for ODEs Derived by the Method of Lines. Euro-Par 2003: 830-839.

[7] Rauber, T., Rünger, G.: Execution Schemes for Parallel Adams Methods. Euro-Par 2004: 708-717.