the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 09 Dec 2025
Automatic tuning of iterative pseudo-transient solvers for modelling the deformation of heterogeneous media
Abstract. Geodynamic modeling has become a crucial tool for investigating the dynamics of Earth deformation across various scales. Such simulations often involve solving mechanical problems with significant material heterogeneities (e.g., strong viscosity contrasts) under nearly incompressible conditions. Recent advancements have enabled the development of iterative solvers based on Dynamic Relaxation or Pseudo-Transient schemes, which require minimal global communication and exhibit quasi-linear scaling on GPU and supercomputing architectures. These solvers incorporate automatic tuning of iterative parameters, including pseudo-time steps and damping coefficients, based on spectral estimates of the discrete operators, ensuring both robust and rapid convergence. We demonstrate the effectiveness of this approach on problems discretized using finite-difference and face-centered finite volume methods, including heterogeneous incompressible Stokes flows. Moreover, the relative algorithmic simplicity of DR-based methods allows for straightforward extensions to compressible flow, multiphase flow, and nonlinear constitutive laws, opening promising avenues for large-scale, high-resolution simulations of geoscientific problems.
Competing interests: one of the authors is also editor
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.- Preprint
(5534 KB) - Metadata XML
- BibTeX
- EndNote
Status: open (until 03 Feb 2026)
- RC1: 'Comment on egusphere-2025-5641', Alexander Minakov, 27 Dec 2025 reply
-
RC2: 'Comment on egusphere-2025-5641', Anthony Jourdon, 12 Jan 2026
reply
I reviewed the manuscript entitled “Automatic tuning of iterative pseudo-transient solvers for modelling the deformation of heterogeneous media” by Duretz et al. The article presents a dynamic relaxation method to automatically compute the parameters (pseudo-time step and damping term) required by the pseudo-transient solvers. First, they present the principle on a simple 1D Poisson case before applying it to a more complex problem involving saddle-point systems such as the Stokes equation, widely used in computational geodynamics. They show the robustness of their approach with community used benchmarks. The method is showed to be performant and make use of GPU hardware which is very state of the art.
The article is very clear, the methods are well explained and motivated, the reader is guided across the entire procedure and points that could rise questions such as preconditioners are addressed and discussed. In my opinion this article is very relevant for the geodynamic/geoscientific modelling community and should be published by Geoscientific Model Development with very minor corrections that could be handled in a proof correction step.
Minor points:Eq. 7: is sub index “i” next to \hat v supposed to be bold?
What is the sub index “i” representing in this Eq.? In Eqs. 8, 9, 10 it is given as the facet index but not in Eq. 7.
Is n_elem the total number of elements in the domain or the number of elements connected to a face i ?
In Eqs. 8, 9, 10, (and likely 7) is n_fac the total number of facets over the domain \Pi or the number of facets per element?Eq. 10 and Eq. 25: I cannot find a definition of \mathbf b. I think I understand this is the source term, but it was named \mathbf f in the strong form in Eq. 1. If this is it, maybe you should reconsider renaming it consistently across the manuscript, if it is not I think it misses a definition.
l. 95, 99: “hybrid velocity vector”, what is special about this vector to be called hybrid?
Eq. 15: You can consider putting \Delta{\bar t} after the fraction as the bar on the t maybe interpreted as a minus sign in front of \Delta u by a fast reader, but this is very minor.Eq. 25: P is a capital letter in the argument of the functional but small caps in the second integral.
Eq. 26 and Alg. 1 (a): in Eq. 26 \mathbf f is used while in Alg. 1 (a) \mathbf b is used. I think you should consider harmonizing the notation to enhance the readability.
Figure 3: I think that if you add an axis x-y-z it will improve the readability of the figure and make it clearer with respect to the use of x, y and z directions in the text.
l. 345: tau_II is written in letters
l. 367: parenthesis not closed
l. 146, 336, 370: How is the number of 100 iterations chosen? Did you experiment with this number?
l. 388-389: strange citation formatting
Citation: https://doi.org/10.5194/egusphere-2025-5641-RC2 -
RC3: 'Comment on egusphere-2025-5641', Anonymous Referee #3, 22 Jan 2026
reply
The manuscript presents a novel, GPU-ready, iterative solver of Powell-Hestenes with inner Dynamic Relaxation velocity solve (PH/DR). As such, it is scientifically significant and is a good fit for GMD. During the review I did not find any major methodological issues. The manuscript is generally well written, and well illustrated. I found some minor shortcomings. I think the examples chosen don't illustrate well the strength of this new method and I think the method description would benefit from being slightly more didactic. The manuscript is in a very good shape, even without addressing these minor shortcomings, therefore I recommend to accept it subject to minor revisions. I don't see the need to review the revised version, approval by the Editor should be sufficient.
Minor shortcoming #1:
In the introduction the novel PH/DR is being compared to the single-loop Pseudo-Transient (PT) solvers. The authors correctly point out that PT solvers struggle with sharp contrasts in material properties. However, I don't see how the new solver performs better in this regard. For the incompressible Stokes inclusion examples, the authors used the same smoothed viscosity field and got visually identical results, including the spurious oscillations at the interface. The other point the authors make is that PT solvers need optimally tuned iteration parameters, which might change depending on the material property fields. This can indeed be problematic for simulations with significant changes in the material property fields over time. However, the majority of the examples in the manuscript are static, where this not an issue. There is a single dynamic example, the visco-elasto-viscoplastic shear banding problem, but for this case no PT based benchmark is provided. I think it would be nice to find and include at least one example where the PT solver struggles and the PH/DR method prevails. This is a suggestion, not a requirement.
Minor shortcoming #2:
I think a more didactic methodological description, e.g. highlighting how the DR method differs from the accelerated PT method, could be beneficial. Maybe this issue will be solved by uploading example codes of the same problem with PT and PH/DR solvers.
Specific comments:
- Line 43: “exact derivation of optimal parameters based on simplified model configuration (Räss et al., 2022; Alkhimenkov et al., 2021)” Add either “a simplified” or “configurations”
- Please always indicate when the provided references are examples by adding e.g. For example, near line 80 "Geometrically, FD has been successfully implemented in Cartesian and polar coordinates (Räss et al., 2017), spherical configurations (Tackley, 2008; Macherel et al., 2024), and on adaptively refined meshes (Gerya, 2013; Goldade et al., 2019)."
- eq. 18 and eq. 19.: c_CFL and c_damp appear without introduction. I guess c_CFL is supposed to be non-negative i.e. 0<c_CFL<1. If this is not the case please specify and elaborate.
- Near line 135: "In either case, two iteration parameters need to be determined as a function of the minimum (λmin) and maximum (λmax) eigenvalues of the discrete problem (Papadrakakis, 1981; Oakley and Knight, 1995):" Please specify which two parameters are you referring to. You give a formula for delta pseudo time for c, but both of them contain an undetermined constant. Please specify how to choose optimal values for c_CFL and c_damp. Also, are these values uniform or are they variables?
- Section 5: Please clearly define all newly introduced variables, such as gamma or r_SCR.
- Algorithm 1: please make it unmistakable that k is the index value for the repeat loop and not the coefficient from eq. 24.
- Algorithm 1: How are a and b are supposed to be initialized? For the first time step no old values of variables and residuals are available.
- Section 6.1., please define the relation between eta_ref, eta_min and eta_max.
Citation: https://doi.org/10.5194/egusphere-2025-5641-RC3
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 197 | 114 | 19 | 330 | 19 | 12 |
- HTML: 197
- PDF: 114
- XML: 19
- Total: 330
- BibTeX: 19
- EndNote: 12
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
In their paper, T. Duretz and colleagues suggest improvements to existing strategies of numerical geodynamic modeling through exploring the application of the direct relaxation (DR) method to various geodynamic problems. The key contributions compared to previous works are in the introduction of automated iteration parameter selection in the pseudo-transient method (pseudo-time stepping and damping) and solving incompressible Stokes flow equations using this method combined with Powell-Hestenes iterations. The iteration parameters are determined based on the eigenvalues of the discrete problem. The paper also addresses challenges associated with large viscosity contrasts and the enforcement of incompressibility, showing a systematic analysis. Finite Difference (FD) and Face-Centered Finite Volume (FCFV) discretization methods are employed, with the FCFV approach yielding smooth solutions across viscosity discontinuities. The study is complemented by numerical examples and an accompanying code repository, providing useful resources for further development of practical applications within and beyond the geodynamic community. The paper is well written, but it can still be improved in a few places with more explanation.
Comments: