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: final response (author comments only)
-
RC1: 'Comment on egusphere-2025-5641', Alexander Minakov, 27 Dec 2025
-
AC1: 'Reply on RC1', Thibault Duretz, 08 Mar 2026
We thank the reviewer for providing his review and suggestions for improving the mansucript. Please find a below a point-by-point answer to the reviewer's comments.
With our best greetings,
The authors###############################
Comments:
1: explain all used notations
–> We have taken care to define all symbols in equation 17: “i” subscript of vhat should be italic
–> This is now correctedLine 95: explain “hybrid” velocity
–> We added a sentence that explains what hybrid velocity means: "The hybrid velocity is an independent variable defined on element faces, different to the element-wise velocity. It serves as a Lagrange multiplier–type trace variable that enforces weak continuity of the velocity field across elements. In HDG and FCFV, the hybrid velocity allows static condensation of the interior degrees of freedom."Line 155: “in” Fig. 2
–> done3: Colormaps are difficult to see
–> we have improved the figure. For some reason all .png do not appear sharp in the rendered pdf. This should not be the case once the document is edited. The color bars and their indications should now appear bigger, in addition the figure resolution should now be appropriate.8: Explain both solid and dotted lines in b) and c). Indicate that x and y are [n.d.].
–> The explanations and indications were added to the caption.9: Can you explain the behavior of fluid continuity? Indicate that x and y are [n.d.].
–> Yes, the initial behaviour is controlled by the initial values of the fluid pressure field. In our case, it's simply 0 and the initial residual is close to 0. Thus, in the first iterations, the pressure residual can only grow. After that, the fluid pressure is updated together with velocity (inner DR loop), it thus exhibits the same sawtooth pattern as velocity. It bumps up every time the total pressure is updated, and then drops again to a lower level, until the next total pressure update. Information about the x and y axis were added to the caption.Line 315: An example with a porosity and permeability range representative of the sedimentary basins can be useful.
–> Yes, we agree. Here, we preferred to keep these model examples in dimensionless units and simple configuration, such that they can be easily reproduced. Application to practical hydromechanical problems will performed in future studies.Line 365: close parentheses after A3
–> doneA1: Indicate that x, y, and P are [n.d.].
–> This is added to the caption.A2: Indicate that x, y, and P are [n.d.]. Why does pressure converge more slowly than velocity in this study?
–> The information was added to the caption. In the PH/DR scheme, the pressure is not updated in the same loop as the velocity. This is a main difference to the approach of Raess et al., 2022. For this reason, the convergence of the pressure residual does not closely follow that of the momentum residuals. Text was added to the appendix.Reference Rass et al. 2017 appears twice in the list
–> thanks for observing this, the double was removed.Citation: https://doi.org/10.5194/egusphere-2025-5641-AC1
-
AC1: 'Reply on RC1', Thibault Duretz, 08 Mar 2026
-
RC2: 'Comment on egusphere-2025-5641', Anthony Jourdon, 12 Jan 2026
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 -
AC2: 'Reply on RC2', Thibault Duretz, 08 Mar 2026
We wish to thank the reviewer for spending time on reviewing our work and for his positive comments towards our manscuript. Below, we include point-by-point answers to the minor points raised by the reviewer.
Best regards,
The authors###############################
Minor points:
Eq. 7: is sub index “i” next to \hat v supposed to be bold?
–> This is now fixedWhat 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.
–> Yes, i is a facet index. In this equation, the index could simply be dropped as the residual is anyway defined at the facet, thus only involving a summation over the 2 connected elements.Is n_elem the total number of elements in the domain or the number of elements connected to a face i ?
–> Here, n_elem is the number of elments connected to the face, which is anyways 2, independently of the number of dimensions. We have now revised the notation.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?
–> n_fac indeed corresponds to the number of facets per element. This precision was added to the text.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.
–> thanks for noticing this, indeed it should be b everywhere to be consistent with the strong form.l. 95, 99: “hybrid velocity vector”, what is special about this vector to be called hybrid?
–> This terminology is used in FCFV context as these velocity dofs are specifically added to enforce velocity continuity across the elements. A clarification was added to the main text.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.
–> Oh, I did not think about that eventuality. We have added a space in between. Now it should be clear that the bar is laying on top of t. Thanks for noting this.Eq. 25: P is a capital letter in the argument of the functional but small caps in the second integral.
–> Thanks for noticing this, it's now fixedEq. 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.
–> this is now harmonized consistentlyFigure 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.
–> yes, that would be a nice addition but I did not manage to generate a nice looking high resolution png export with a triad using paraview.l. 345: tau_II is written in letters
–> This is now fixedl. 367: parenthesis not
–> Fixedl. 146, 336, 370: How is the number of 100 iterations chosen? Did you experiment with this number?
–> yes, we did experiment and we empirically found this number to be a good compromise. The aim is not to do these operations (reductions) too frequently. From our experience with different benchmarks and geodynamic miniapps, n is typically 10<n<400. We added a sentence around line 172 to clarify that point.l. 388-389: strange citation formatting
fixedCitation: https://doi.org/10.5194/egusphere-2025-5641-AC2
-
AC2: 'Reply on RC2', Thibault Duretz, 08 Mar 2026
-
RC3: 'Comment on egusphere-2025-5641', Anonymous Referee #3, 22 Jan 2026
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 -
AC3: 'Reply on RC3', Thibault Duretz, 08 Mar 2026
We thank the reviewer for the careful evaluation of our manuscript. The constructive suggestions, remarks, and corrections have been addressed in the revised version. We believe these revisions have strengthened the manuscript and hope that they satisfactorily address the reviewer’s concerns.
Best regards,
The authors###############################
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.
–> The PH/DR solver does not address issues related to the spatial discretization. In the presence of sharp material contrasts, and without a conforming mesh or an interface upscaling technique, spurious oscillations at material interfaces are to be expected. In terms of the total number of iterations and scaling behavior, the proposed solver performs better than a pseudo-transient (PT) solver that is optimally tuned. However, the primary strength of the method lies in its automatic tuning capability. In contrast to PT approaches, the damping and time-stepping parameters do not need to be derived analytically (when feasible) or adjusted manually. This is particularly advantageous for dynamically evolving problems, where the spectral properties of the operator change over time. Such variations pose a significant challenge for PT-based methods, as the parameters may leave the region of optimality, leading to stagnation or even divergence. The proposed auto-tuning strategy tracks these changes in spectral properties and adjusts the iterative parameters accordingly. Therefore, the key benefit of the method is improved robustness and flexibility, rather than a strict reduction in iteration count relative to an optimally tuned PT solver. We agree that adding a set of practical and more dynamic examples is a valuable suggestion. In the revised manuscript, we included a comparison with models performed using JustRelax. These examples are relevant, as they involve time-dependent geodynamic simulations. All simulations were carried out within the same geodynamic modeling framework (JustRelax.jl), which implements both a static pseudo-transient solver and the PH/DR algorithm described in this work. In these cases, the performance gain achieved by using PH/DR is substantial.
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.
–> We agree with the reviewer and have added a specific appendix section that explains, in descriptive terms, the main similarities and differences between the PT and the DR and PHDR solution strategies. In addition, Appendix 4 provides a quantitative comparison of the two approaches, based on the reference model presented in Raess et al. (2022).
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”
–> thanks, the "s" was indeed missing -
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)."
-> we added the missing e.g. -
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.
-> we do now introduce and elaborate on c_CFL and c_damp. -
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?
-> we have added more details about the parameter selection procedure. The real challenge is the estimation of the extremal eigenvalues. c_CFL and c_damp are single scalar parameters, they do not vary over space. c_CFL should be maximised (0.99) for optimal convergence, while the other c_damp can be varied. In practice there is no need to vary it much from the predicted optimal value of 1 (Oakley, 1995). -
Section 5: Please clearly define all newly introduced variables, such as gamma or r_SCR.
-> gamma is introduced line 184. r_scr is now introduced as "the residual of the velocity Schur complement" -
Algorithm 1: please make it unmistakable that k is the index value for the repeat loop and not the coefficient from eq. 24.
-> thanks for pointing this out. This is now fixed. -
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.
-> This is now clarified in the section which presents the estimation of lambda_min. -
Section 6.1., please define the relation between eta_ref, eta_min and eta_max.
-> This is now explicitely stated
Citation: https://doi.org/10.5194/egusphere-2025-5641-AC3 -
-
AC3: 'Reply on RC3', Thibault Duretz, 08 Mar 2026
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 376 | 220 | 33 | 629 | 25 | 17 |
- HTML: 376
- PDF: 220
- XML: 33
- Total: 629
- BibTeX: 25
- EndNote: 17
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: