A log-linear time algorithm for the elastodynamic boundary integral equation method

This paper introduces a fast and memory-efficient algorithm called FDP=H-matrices that utilizes fast domain partitioning and plane-wave approximations to solve transient elastodynamic boundary integral equations with $\mathcal{O}(N \log N)$ memory and $\mathcal{O}(NM \log N)$ computation time, significantly reducing the costs of traditional time-marching implementations while maintaining accuracy.

The Gist

Imagine you are trying to predict how a massive earthquake will shake a city. To do this, scientists use a mathematical tool called the Boundary Integral Equation Method (BIEM). Think of this method as a way to simulate the earthquake by only looking at the "skin" of the problem (the fault lines and the ground surface) rather than the entire 3D volume of the Earth. This is efficient because it ignores the empty space inside the ground.

However, there's a huge catch. In the traditional way of doing this, the computer has to remember every single interaction between every point on the fault and every point on the ground, for every single moment in time.

The Problem: The "Infinite Library"

Imagine you have a library where you need to find a specific book.

• The Old Way: You have to walk to every single shelf (every point on the fault), check every single book (every time step), and write down a note about how that book affects your current location.
• The Cost: If you have 1,000 points and 1,000 time steps, you have to write down 1,000,000 notes just for one moment. If you want to simulate a whole earthquake, the number of notes explodes to trillions.
• The Result: Your computer runs out of memory (RAM) and takes years to finish the calculation. It's like trying to carry the entire Library of Congress in your backpack.

The Solution: FDP=H-Matrices

The authors of this paper developed a new algorithm called FDP=H-matrices. Think of this as a "Smart Librarian" that doesn't just carry books; it understands the story of the earthquake so well that it can summarize thousands of notes into a single sentence.

Here is how they did it, broken down into four clever tricks:

1. The "Wavefront" Partition (FDPM)

Earthquakes send out waves (like ripples in a pond). These waves travel at specific speeds.

• The Trick: The new method divides time into three zones:
  1. The Wave Zone: The exact moment the wave hits.
  2. The Aftermath: The time between the first wave and the second.
  3. The Calm: The time after everything has settled.
• Why it helps: In the "Aftermath" and "Calm" zones, the physics becomes simple and predictable. The complex math can be simplified into a "Space part" (where things are) and a "Time part" (when things happen), which are easy to separate.

2. The "Compression" (H-Matrices)

In the "Wave Zone," things are messy and singular (like a sudden shock). Usually, computers can't compress this data.

• The Trick: The authors realized that even though the wave is sharp, the shape of the wave as it travels follows a predictable pattern (it gets weaker as it spreads out, like a flashlight beam).
• The Analogy: Imagine a crowd of people shouting. Instead of recording every individual voice, you record the "average volume" and the "direction" of the shout. You can describe the whole crowd's noise with just two numbers.
• The Result: They used a technique called H-matrices to compress the messy wave data into a tiny, low-rank summary. Instead of storing millions of numbers, they store a few hundred.

3. The "Plane Wave" Shortcut (ART)

When waves travel far away, they look flat (like a plane wave).

• The Trick: The algorithm assumes that for distant groups of points, the wave arrives at almost the same time. It uses a "relay station" (a representative point) to calculate the travel time for the whole group.
• The Analogy: Instead of calculating the exact time it takes for a message to reach every person in a stadium individually, you calculate the time to reach the center of the stadium and add a tiny correction for people sitting in the front vs. the back.
• The Result: This removes the need to store the history of every single point, saving massive amounts of memory.

4. The "Staircase" Sampling (Quantization)

In the "Aftermath" zone, the waves change slowly over time.

• The Trick: Instead of checking the wave every millisecond, the algorithm checks it every second, then every 10 seconds, then every 100 seconds. It takes "steps" up the time ladder.
• The Analogy: If you are watching a slow-moving glacier, you don't need a photo every second. You can take a photo every hour, and it will still look accurate.
• The Result: This drastically cuts down the number of calculations needed for the long tail of the earthquake simulation.

The Grand Result

By combining these four tricks, the authors turned a problem that was impossible for standard computers into one that is fast and cheap.

• Old Method: Memory usage grows like $N^2$ (if you double the points, you need 4x the memory). Time grows like $N^2$ (if you double the points, it takes 4x longer).
• New Method: Memory usage grows like $N \log N$ (if you double the points, you need barely a little more memory). Time grows similarly.

In simple terms:
If the old method was like trying to carry a mountain of sand in your hands, the new method is like using a vacuum cleaner that sucks up the sand and compresses it into a tiny, lightweight bag. You can now simulate massive, complex earthquakes on a standard laptop that used to require a supercomputer, and you can do it in minutes instead of years.

This breakthrough allows scientists to model disasters with much higher detail and accuracy, potentially helping us build safer cities and better understand the Earth's dynamics.

Technical Summary

1. Problem Statement

The paper addresses the computational bottleneck associated with the Spatiotemporal Boundary Integral Equation Method (ST-BIEM) for transient elastodynamic problems (e.g., dynamic rupture, wave scattering).

• The Challenge: In traditional ST-BIEM, the solution requires convolving boundary variables (slip/opening rates) with a dense kernel tensor over $N$ boundary elements and $M$ time steps.
• Complexity: The naive implementation results in:
  • Memory: $O(N^2 M)$ to store the discrete kernel and $O(NM)$ for boundary history.
  • Computation Time: $O(N^2 M^2)$ for a single run (or $O(N^2 M)$ per time step).
• Limitations of Existing Methods:
  • Volume-based methods (FDM/FEM): Scale as $O(N_v M)$ but require discretizing the entire volume ($N_v \gg N$), losing the efficiency of boundary-only discretization.
  • Plane-Wave Time-Domain (PWTD): Reduces time to $O(NM \log^2 N)$ but requires $O(NM)$ memory to store boundary histories and involves complex analytical formulations.
  • Hierarchical Matrices (H-matrices): Efficient for static problems ($O(N \log N)$) but struggle with elastodynamic kernels because the singularities (wavefronts) are distributed along the wave arrival time in the space-time domain, preventing standard low-rank approximations.
  • Convolution Quadrature Methods (CQM): Efficient in the frequency domain but often rely on high-frequency approximations that struggle with low-frequency motions.

2. Methodology: FDP=H-Matrices

The authors propose a novel algorithm named Fast Domain Partitioning Hierarchical Matrices (FDP=H-matrices). This method combines four distinct modules to achieve $O(N \log N)$ scaling in both memory and time per step.

A. Core Architecture

The algorithm integrates:

1. Fast Domain Partitioning Method (FDPM):

   • Partitions the space-time domain of the kernel into three regions based on wave arrival times:
     • Domain F: Contains the impulsive P- and S-wavefronts (singularities).
     • Domain I: The region between P- and S-waves (near-field term).
     • Domain S: The static equilibrium region after S-waves.
   • In Domains I and S, the kernel naturally separates into spatial and temporal components, allowing for low-rank approximation.

2. H-Matrices (Hierarchical Matrices):

   • Applied to the spatial dependence of the kernels in all domains.
   • Innovation: The authors apply H-matrices to Domain F (the wavefronts). By integrating the kernel over the wavefront duration (creating an "amplitude term"), they transform the impulsive delta functions into smooth, geometrically spreading functions that admit low-rank approximation (rank $O(1)$).

3. Average Reduced Time (ART):

   • A plane-wave approximation technique used to handle the travel time ($t_{ij}$) dependency in Domain F.
   • Approximates the travel time between a source $j$ and receiver $i$ as $t_{ij} \approx \bar{t}_j + \delta t_i$, separating source and receiver dependencies.
   • This allows the algorithm to avoid storing the full history of boundary variables, a major memory bottleneck in previous methods.

4. Quantization:

   • A sparse resampling technique applied to the temporal integration in Domain I.
   • It approximates the time-dependent part of the kernel using a piecewise-constant staircase function, reducing the number of time steps required for convolution without significant loss of accuracy.

B. Arithmetic and Memory Efficiency

• Memory Reduction: The algorithm eliminates the need to store the $O(NM)$ history of boundary variables. Instead, it uses a recursive "dynamic programming" approach to update a "representative stress" vector.
  • Total Memory: $O(N \log N)$.
• Time Complexity:
  • Per Time Step: $O(N \log N)$.
  • Total Run: $O(NM \log N)$.
• Two Schemes: The authors propose two admissibility conditions for the ART:
  • Constant $\eta$ scheme: Achieves $O(N \log N)$ scaling.
  • Constant $\eta^2 d_{ist}$ scheme: Achieves $O(N^{3/2})$ scaling but offers higher accuracy for distant clusters.

3. Key Contributions

1. First Log-Linear Elastodynamic Solver: This is the first versatile algorithm to achieve $O(N \log N)$ total memory and $O(N \log N)$ time complexity per time step for transient elastodynamic problems on arbitrary geometries.
2. Solving the Singularity Problem: Successfully adapts H-matrices to the elastodynamic wave equation by integrating the kernel over the wavefront (Domain F), effectively handling the distributed singularities that previously hindered H-matrix application.
3. Memory-Efficient Arithmetic: Developed a new arithmetic scheme that removes the $O(NM)$ memory requirement for boundary history, a significant improvement over PWTD and standard time-marching BIEM.
4. Analytical and Numerical Validation: Provided rigorous error bounds for the ART and Quantization modules and validated them through extensive numerical experiments.

4. Results

The authors validated the method using 2D anti-plane dynamic rupture simulations (with extensions to 3D discussed theoretically).

• Cost Reduction:
  • Compared to the original ST-BIEM ($O(N^2 M)$ memory/time), FDP=H-matrices demonstrated a reduction to $O(N \log N)$.
  • Numerical experiments showed a memory and time reduction factor of approximately $N M / \log N$ compared to orthodox implementations.
• Accuracy:
  • Dynamic Rupture: The method reproduced rupture propagation speeds and slip-rate patterns with high fidelity.
  • Error Magnitude: Relative errors were consistently below 0.4% (and often <0.3%) even after hundreds of time steps, which is significantly smaller than typical numerical dispersion errors in standard BIEM.
  • Parameter Sensitivity: The method remained robust across varying parameters ($\epsilon_{ACA}$, $\eta$, $l_{min}$). The constant $\eta$ scheme showed non-dispersive wave speed errors.
• Scalability: The cost scaling held true for both planar and non-planar (kinked) fault geometries.

5. Significance

• Enabling Large-Scale Simulations: By reducing the memory bottleneck from $O(N^2 M)$ to $O(N \log N)$, this algorithm allows researchers to simulate much larger problems (larger $N$ and $M$) on standard hardware, or simulate the same problems with significantly fewer resources.
• Versatility: Unlike frequency-domain methods or specific spectral methods, FDP=H-matrices work for arbitrary boundary geometries and general transient elastodynamic problems.
• Future Applications: The method opens the door for high-fidelity simulations of complex earthquake rupture dynamics, seismic wave scattering in complex geological structures, and other wave phenomena where boundary integral methods are preferred over volume methods.
• Theoretical Advancement: It bridges the gap between the efficiency of H-matrices and the physics of wave propagation, offering a new paradigm for solving hyperbolic partial differential equations via boundary integrals.

In conclusion, the paper presents a breakthrough in computational seismology and elastodynamics, transforming the ST-BIEM from a computationally prohibitive method for large-scale transient problems into a highly efficient, scalable, and accurate tool.