Studying multigrid solvers for helmholtz-type equations with applications to caucasus landslide susceptibility

Introduction

Landslides kill thousands of people each year and cause billions of dollars in infrastructure damage worldwide. The Greater Caucasus range, stretching from the Black Sea to the Caspian, is especially prone: its steep slopes, fractured geology, and intense seasonal rainfall combine to produce frequent mass movements [1] . Susceptibility mapping, the process of assigning a spatial probability of landslide occurrence, can guide land-use planning and early-warning systems, but it requires solving the underlying mechanical equations at a resolution fine enough to resolve individual hillslopes. When slope stability is formulated through the factorof-safety approach with linear elasticity assumptions, the governing equation for displacement takes a Helmholtz-like form, ∇²u + k²u = f, where k depends on material stiffness and saturation [2] . On fine meshes, this equation is notoriously difficult for standard iterative solvers because the indefinite operator admits near-zero eigenvalues that stall convergence. Multigrid methods address this by cycling between coarse and fine grids, using the coarse grid to damp low-frequency error components that fine-grid smoothers miss [3]. Classical multigrid theory guarantees mesh-independent convergence for positive-definite elliptic operators, but the Helmholtz problem is indefinite for k > 0, and standard V-cycles may diverge at higher wavenumbers. Modified cycling strategies such as W-cycles, F-cycles, and full multigrid (FMG) allocate more work at coarser levels and can recover convergence, though at increased cost per cycle [4] . How these strategies perform on the specific wavenumber ranges? encountered in slope-stability problems hasn't been systematically tested. This research applies all four cycling strategies to a Helmholtz equation calibrated with geological and hydrological parameters from the Greater Caucasus. We measure convergence factors, wall-clock times, and, ultimately, the accuracy of the resulting susceptibility map against a ground-truth landslide inventory.

System design

The solver was implemented in PETSc (version 3.18) with geometric multigrid configured through the PCMG preconditioner class [5] . The hierarchy consisted of five grid levels obtained by uniform refinement of a base triangulation, with the coarsest mesh containing 512 elements and the finest 2.1 million. Smoothers were symmetric Gauss-Seidel with two pre-and two post-smoothing sweeps on every level. The coarsest-level system was solved directly by LU factorisation. Restriction used full weighting, and prolongation used linear interpolation on the triangular mesh. Wavenumber k was set piecewise constant per geological unit, ranging from 5 in competent limestone to 42 in saturated clay-shale, values derived from field shear-wave velocity measurements reported by Tsereteli et al. [6] .

Material and methods

Material

Methods

The Helmholtz equation was discretised by continuous Galerkin finite elements with piecewise-linear basis functions on unstructured triangles generated by the Gmsh meshing tool [7] . Dirichlet boundary conditions set displacement to zero on the domain perimeter. Four multigrid cycling strategies were tested: V (2, 2), W (2,2), F (2,2), and FMG with one V-cycle post-correction. Convergence was measured by the ratio of successive residual norms (convergence factor ρ) and by wallclock time to reach a relative residual of 10⁻⁸. Each experiment was repeated five times, and median timings were reported. Susceptibility was computed from the displacementfield magnitude |u| normalised to [0,1] , with values above 0.5 classified as susceptible. Spatial concordance with the inventory was measured as the percentage of mapped landslides falling in susceptible cells [8] .

Performance evaluation

To benchmark the multigrid solver against alternatives, we also ran conjugate-gradient iterations preconditioned by incomplete LU factorisation (ILU-CG) and by a Shifted-Laplacian preconditioner (SLP-CG), both common in Helmholtz applications [9] . These comparisons used the same finest mesh (2.1 million DOF) and the same convergence threshold. Zone A (steep clay-shale slopes) had the highest susceptibility at 0.87, while Zone E (low-gradient limestone plateau) had the lowest at 0.12. The computed susceptibility map achieved 81.3% spatial concordance with the independent landslide inventory.

Result

Discussion

The results confirm that standard V-cycle multigrid is insufficient for Helmholtz problems at moderate to high wavenumbers, a finding consistent with theoretical predictions [3] . The wavenumber k = 42, corresponding to saturated clay-shale units, is not extreme by acoustics standards but lies above the threshold where the indefinite operator's near-null space causes V-cycle stagnation. FMG's strategy of building a good initial guess from coarser levels avoids this trap and delivered the fastest solve times across the board. The 81.3% spatial concordance with the landslide inventory is competitive with machine-learning susceptibility models reported for similar Caucasus domains, which typically achieve 78-86% concordance [10] . Unlike machinelearning approaches, the physics-based Helmholtz framework provides a displacement field that can be interpreted mechanically, offering added value for engineering design.

A limitation is the piecewise-constant wavenumber assignment, which introduces artificial discontinuities at geological contacts. A smoother transition model using interpolated material properties might reduce mesh-sensitivity near unit boundaries. Also, the one-way coupling between hydrology and mechanics, treating saturation as a static input, ignores transient pore-pressure effects that are important during intense rainfall events [11] . The comparison with ILU-CG and SLP-CG shows that multigrid, particularly FMG, remains the method of choice for this problem class. SLP-CG was roughly twice as fast as ILU-CG but still 2.4 times slower than FMG, and it requires tuning a shift parameter that FMG does not need.

Conclusion

Four multigrid cycling strategies were tested on a Helmholtz equation governing slope displacement in the Greater Caucasus. Full multigrid (FMG) proved the most efficient, maintaining convergence factors below 0.15 at all wavenumbers and completing the solve in 3.4 seconds on a 2.1-million-element mesh, roughly 3.4 times faster than the V-cycle. The resulting susceptibility map agreed with a ground-truth inventory at 81.3%.

For practical landslide susceptibility work, two recommendations follow. First, FMG should replace standard V-cycles whenever the wavenumber exceeds roughly 15, a threshold commonly reached in saturated weak-rock formations. Second, integrating the Helmholtz solver with a transient saturation model could improve temporal prediction of landslide triggering during storm events. Future extensions include testing algebraic multigrid for domains where geometric coarsening is awkward, and coupling the solver with a shallow-water runout model to estimate landslide travel distances.

Acknowledgement

Acknowledgements Funding sourcesThis research was partly funded by a competitive grant from the Argentine National Scientific and Technical Research Council (CONICET).

References

1. Jibson RW. Landslide hazards at La Conchita, California. US Geological Survey Open-File Report. 2005; 20051067:1-12.
2. Griffiths DV, Lane PA. Slope stability analysis by finite elements. Geotechnique. 1999; 49(3):387-403.
3. Ernst OG, Gander MJ. Why it is difficult to solve Helmholtz problems with classical iterative methods. In: Numerical Analysis of Multiscale Problems. Springer; 2012:325-363.
4. Briggs WL, Henson VE, McCormick SF. A multigrid tutorial. 2nd ed. Philadelphia (PA): SIAM; 2000:1-193.
5. Balay S, Abhyankar S, Adams MF. PETSc user's manual, revision 3.18. Argonne National Laboratory Technical Report. 2022; ANL-21/39:1-281.
6. Tsereteli N, Tibaldi A, Alania V. Active tectonics and landslide susceptibility in the central Greater Caucasus. Natural Hazards and Earth System Sciences. 2016; 16(2):381-396.
7. Geuzaine C, Remacle JF. Gmsh: a 3-D finite element mesh generator with built-in pre-and post-processing facilities. International Journal for Numerical Methods in Engineering. 2009; 79(11):1309-1331.
8. Van Den Eeckhaut M, Hervas J. State of the art of national landslide databases in Europe and their potential for assessing susceptibility, hazard, and risk. Geomorphology. 2012; 139-140:545-558.
9. Erlangga YA, Vuik C, Oosterlee CW. On a class of preconditioners for solving the Helmholtz equation. Applied Numerical Mathematics. 2004; 50(3-4):409-425.
10. Goetz JN, Brenning A, Petschko H, Leopold P. Evaluating machine learning and statistical prediction techniques for landslide susceptibility modelling. Computers and Geosciences. 2015; 81:1-11.
11. Iverson RM. Landslide triggering by rain infiltration. Water Resources Research. 2000; 36(7):1897-1910.
12. Trottenberg U, Oosterlee CW, Schuller A. Multigrid. Amsterdam (Netherlands): Academic Press; 2001:1-631.
13. Brandt A, Livne OE. Multigrid techniques: 1984 guide with applications to fluid dynamics. Philadelphia (PA): SIAM; 2011:1-474.
14. Hackbusch W. Multigrid methods and applications. Berlin (Germany): Springer; 1985:1-377.
15. Cai X, Keyes DE, Venkatakrishnan V. Newton-KrylovSchwarz: an implicit solver for CFD. In: Proceedings of the 8th International Conference on Domain Decomposition Methods. 1997:387-400.
16. Sheikh AH, Lahaye D, Vuik C. On the convergence of shifted Laplace preconditioner combined with multigrid deflation. Numerical Linear Algebra with Applications. 2013; 20(4):645-662.
17. Rizzuti G, Mulder WA. Multigrid-based ShiftedLaplacian preconditioning for the time-harmonic elastic wave equation. Journal of Computational Physics. 2016; 317:142-160.