the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Benchmarking the accuracy of higher order particle methods in geodynamic models of transient flow
Abstract. Numerical models are a powerful tool for investigating the dynamic processes in the interior of the Earth and other planets, but the reliability and predictive power of these discretized models depends on the numerical method, as well as an accurate representation of material properties in space and time. In the specific context of geodynamic models, particle methods have been applied extensively because of their suitability for advection-dominated processes, and have been used in applications such as tracking the composition of solid rock and melt in the Earth’s mantle, fluids in lithospheric- and crustal- scale models, light elements in the liquid core, and deformation properties like accumulated finite strain or mineral grain size, along with many applications outside the Earth sciences.
There have been significant benchmarking efforts to measure the accuracy and convergence behavior of particle methods, but these efforts have largely been limited to instantaneous solutions, or time-dependent models without analytical solutions. As a consequence, there is little understanding about the interplay of particle advection errors and errors introduced in the solution of the underlying transient, nonlinear flow equations. To address these limitations, we present two new dynamic benchmarks for transient Stokes flow with analytical solutions that allow us to quantify the accuracy of various advection methods in nonlinear flow. We use these benchmarks to measure the accuracy of our particle algorithm as implemented in the ASPECT geodynamic modeling software against commonly employed field methods and analytical solutions. In particular, we quantify if an algorithm that is higher-order accurate in time will allow for better overall model accuracy and verify that our algorithm reaches its intended optimal convergence rate. We then document that the observed increased accuracy of higher-order algorithms matters for geodynamic applications with an example of modeling small-scale convection underneath an oceanic plate and show that the predicted place and time of onset of small-scale convection depends significantly on the chosen particle advection method.
Descriptions and implementations of our benchmarks are openly available and can be used to verify other advection algorithms. The availability of accurate, scalable and efficient particle methods as part of the widely used open source code ASPECT will allow geodynamicists to accurately investigate more complex time-dependent geodynamic processes, such as elastic deformation, anisotropic fabric development, melt generation and migration, and grain damage.
-
Notice on discussion status
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
-
Preprint
(5240 KB)
-
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(5240 KB) - Metadata XML
- BibTeX
- EndNote
- Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2023-2765', Anonymous Referee #1, 05 Jan 2024
This is a well-written technical paper that analyses the numerical errors of different advection methods commonly used in geodynamic modelling and compares them to higher-order particle methods. The paper also presents two very nice new advection benchmarks for transient flow, which highlight the importance of not only benchmarking but also developing new benchmarks. A pdf with some minor comments and suggestions is attached.
-
RC2: 'Comment on egusphere-2023-2765', Anonymous Referee #2, 09 Jan 2024
General comments
I have read the manuscript of Gassmöller et al. with interest. The study investigates the accuracy of the particle in celle methods, widely used in geodynamic modelling.
New analytical solutions are provided (and used) and I'm sure they will be very useful for the community.
The main conclusion is that geodynamic simulators should include second-order in time particle advection solvers.
I think the manuscript is worth publication after addressing the moderate comments listed below.
Addressing these comments should help with:
- providing a better overview of existing literature
- clarifying whether their conclusions are applicable in the general case (Main comment #1)
- demonstrating their points in a more convincing manner (Main comment #2)
section 1
Reading the first part of the introduction gives the feeling that geodynamic modelling started in the lated 2010s. I would suggest to also cite older work. Numerical geodynamic models can be traced back to the late 1970s.
l 36. you may want to oppose Eulerian to Lagrangian schemes (instead of field vs particle-based schemes).
l 37-46 the intro would better flow if this paragraph would be better located elsewhere (e.g. section 3). As such it reads weird. One would expect the intro to focus on the topic of advection (which is introduce in the line above this paragraph and discussed right after)
l 58 yes, this is right, there is no systematic study of this type of error. Some examples are shown in Gerya and Yuen 2003 (e.g. Fig 12) and a potential solution is proposed (subgrid diffusion), the scheme is further extended to visco-elasticity in Gerya and Yuen (2007). This could be mentioned.
I would mention Weinberg & Schmeling (1992) that uses particles with RK4 and a predictor-corrector scheme
https://www.sciencedirect.com/science/article/abs/pii/0191814192901034
section 2
eq 3 I’m not sure there is a point to include a diffusion term which is not used later on. Otherwise inactive terms could also appear in other governing equations (inertia, compressibility…).
l 98 shear heating is not restricted to compressible flow. You probably want to refer to adiabatic heating instead.
section 3
I would definitely put the paragraph about Aspect history here (it’s currently in the middle of the introduction).
In the description of the each scheme (FE, RK2, RK2FOT) I would clearly indicated how many full (non-linear) mechanical solves are needed to integrate one time step.
By reading the last paragraph, it may be worth mentioning somewhere the work of Furuichi and May (2015) on implicit advection:
https://www.sciencedirect.com/science/article/pii/S0010465515000569
There is also intersecting work on time integration of non-linear Stokes problems in Adamuszek et al (2016):
https://www.sciencedirect.com/science/article/abs/pii/S0191814116300013
section 4.1
I think this both the analytic solution and its description are very useful.
section 4.2
Here it reads a bit weird in the sense that the authors describe a modification of a solution used in a previous study and the difficulties encountered. As such it reads a bit like personal notes. Maybe you could state the complete solution in Cartesian coordinates directly and simply mention that it’s different than the one used in the previous study? This would decrease the dependance of the current manuscript on this preceding study.
eq 19 is already defined in the text of section 4.1.
section 5.1
l. 308 here it’s not clear why RK2 delivers third order convergence for velocity.
l 333 not sure what the last sentence actually means.
section 5.1
MAIN COMMENT #1:
l. 365 Here you argue that the spatial accuracy of the mechanical solution and the temporal accuracy of the particle advection scheme should be on the same order of accuracy, right?
If yes, then one should opt for RK2 since the temporal resolution matches the spatial accuracy of a Q2 element.
My worry is that, in practice, model configurations generally include material boundaries with coefficient jumps (e.g. viscosity). In these cases, the order of accuracy of the mechanical solutions (L1, L2) will degrade (down to 1st order?) unless the jump is actually meshed. This last option is most often not selected in geodynamic models that rely on particles. Therefore, despite using second order stencils in space, the solutions are effectively of lower order. These cases are the most generally relevant but are not accounted by the analytical solution nor in the next section.
So, in practice, is it still beneficial to use a second order time integration order if spatial accuracy is effectively of lower order?
section 6
MAIN COMMENT #2:
The comparison here is interesting but this section is bit weak to my opinion.
- Is there any random noise on the initial marker position? Is it exactly the same initial condition for both models?
- The RF2FOT is first order but it should also converge to a physical solution if time step is refined. Ideally, the same pattern as for RK2, right?
- What is missing here is a temporal resolution test for the 2 schemes. It would be nice to see wether RF2FOT converges to similar patter as RK2 if time step is refined.
In fact it’s not even clear if the RK2 is close to any physical solution here. The addition of a time convergence study would greatly improve the content of this study.
Moreover, if the 2 schemes converge to the same physical solution, you could clearly evaluate the computing effort needed to reach the “converged” solution pattern with both methods. This could clearly show the benefit of RK2 on a practical case.
Ideally, an application that includes material jump would have been ideal (e.g. D’’, magma body, lithosphere stratification) but I realise it might be out of the scope of the current study.
p 17 Physical units need not to be italicised.
Citation: https://doi.org/10.5194/egusphere-2023-2765-RC2 -
AC1: 'Comment on egusphere-2023-2765', Rene Gassmoeller, 13 Mar 2024
Attached, please find a reply to the points raised by the reviewers. Reviewer comments are in blue, our replies in black. Issues brought up in reviews but not specifically discussed have simply been fixed in the revision as suggested by the reviewers.
We would like to express our gratitude to the reviewers for the very careful reading of our manuscript, along with the thoughtful comments. Peer review makes papers better, and this observation applies to the current one as well.
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2023-2765', Anonymous Referee #1, 05 Jan 2024
This is a well-written technical paper that analyses the numerical errors of different advection methods commonly used in geodynamic modelling and compares them to higher-order particle methods. The paper also presents two very nice new advection benchmarks for transient flow, which highlight the importance of not only benchmarking but also developing new benchmarks. A pdf with some minor comments and suggestions is attached.
-
RC2: 'Comment on egusphere-2023-2765', Anonymous Referee #2, 09 Jan 2024
General comments
I have read the manuscript of Gassmöller et al. with interest. The study investigates the accuracy of the particle in celle methods, widely used in geodynamic modelling.
New analytical solutions are provided (and used) and I'm sure they will be very useful for the community.
The main conclusion is that geodynamic simulators should include second-order in time particle advection solvers.
I think the manuscript is worth publication after addressing the moderate comments listed below.
Addressing these comments should help with:
- providing a better overview of existing literature
- clarifying whether their conclusions are applicable in the general case (Main comment #1)
- demonstrating their points in a more convincing manner (Main comment #2)
section 1
Reading the first part of the introduction gives the feeling that geodynamic modelling started in the lated 2010s. I would suggest to also cite older work. Numerical geodynamic models can be traced back to the late 1970s.
l 36. you may want to oppose Eulerian to Lagrangian schemes (instead of field vs particle-based schemes).
l 37-46 the intro would better flow if this paragraph would be better located elsewhere (e.g. section 3). As such it reads weird. One would expect the intro to focus on the topic of advection (which is introduce in the line above this paragraph and discussed right after)
l 58 yes, this is right, there is no systematic study of this type of error. Some examples are shown in Gerya and Yuen 2003 (e.g. Fig 12) and a potential solution is proposed (subgrid diffusion), the scheme is further extended to visco-elasticity in Gerya and Yuen (2007). This could be mentioned.
I would mention Weinberg & Schmeling (1992) that uses particles with RK4 and a predictor-corrector scheme
https://www.sciencedirect.com/science/article/abs/pii/0191814192901034
section 2
eq 3 I’m not sure there is a point to include a diffusion term which is not used later on. Otherwise inactive terms could also appear in other governing equations (inertia, compressibility…).
l 98 shear heating is not restricted to compressible flow. You probably want to refer to adiabatic heating instead.
section 3
I would definitely put the paragraph about Aspect history here (it’s currently in the middle of the introduction).
In the description of the each scheme (FE, RK2, RK2FOT) I would clearly indicated how many full (non-linear) mechanical solves are needed to integrate one time step.
By reading the last paragraph, it may be worth mentioning somewhere the work of Furuichi and May (2015) on implicit advection:
https://www.sciencedirect.com/science/article/pii/S0010465515000569
There is also intersecting work on time integration of non-linear Stokes problems in Adamuszek et al (2016):
https://www.sciencedirect.com/science/article/abs/pii/S0191814116300013
section 4.1
I think this both the analytic solution and its description are very useful.
section 4.2
Here it reads a bit weird in the sense that the authors describe a modification of a solution used in a previous study and the difficulties encountered. As such it reads a bit like personal notes. Maybe you could state the complete solution in Cartesian coordinates directly and simply mention that it’s different than the one used in the previous study? This would decrease the dependance of the current manuscript on this preceding study.
eq 19 is already defined in the text of section 4.1.
section 5.1
l. 308 here it’s not clear why RK2 delivers third order convergence for velocity.
l 333 not sure what the last sentence actually means.
section 5.1
MAIN COMMENT #1:
l. 365 Here you argue that the spatial accuracy of the mechanical solution and the temporal accuracy of the particle advection scheme should be on the same order of accuracy, right?
If yes, then one should opt for RK2 since the temporal resolution matches the spatial accuracy of a Q2 element.
My worry is that, in practice, model configurations generally include material boundaries with coefficient jumps (e.g. viscosity). In these cases, the order of accuracy of the mechanical solutions (L1, L2) will degrade (down to 1st order?) unless the jump is actually meshed. This last option is most often not selected in geodynamic models that rely on particles. Therefore, despite using second order stencils in space, the solutions are effectively of lower order. These cases are the most generally relevant but are not accounted by the analytical solution nor in the next section.
So, in practice, is it still beneficial to use a second order time integration order if spatial accuracy is effectively of lower order?
section 6
MAIN COMMENT #2:
The comparison here is interesting but this section is bit weak to my opinion.
- Is there any random noise on the initial marker position? Is it exactly the same initial condition for both models?
- The RF2FOT is first order but it should also converge to a physical solution if time step is refined. Ideally, the same pattern as for RK2, right?
- What is missing here is a temporal resolution test for the 2 schemes. It would be nice to see wether RF2FOT converges to similar patter as RK2 if time step is refined.
In fact it’s not even clear if the RK2 is close to any physical solution here. The addition of a time convergence study would greatly improve the content of this study.
Moreover, if the 2 schemes converge to the same physical solution, you could clearly evaluate the computing effort needed to reach the “converged” solution pattern with both methods. This could clearly show the benefit of RK2 on a practical case.
Ideally, an application that includes material jump would have been ideal (e.g. D’’, magma body, lithosphere stratification) but I realise it might be out of the scope of the current study.
p 17 Physical units need not to be italicised.
Citation: https://doi.org/10.5194/egusphere-2023-2765-RC2 -
AC1: 'Comment on egusphere-2023-2765', Rene Gassmoeller, 13 Mar 2024
Attached, please find a reply to the points raised by the reviewers. Reviewer comments are in blue, our replies in black. Issues brought up in reviews but not specifically discussed have simply been fixed in the revision as suggested by the reviewers.
We would like to express our gratitude to the reviewers for the very careful reading of our manuscript, along with the thoughtful comments. Peer review makes papers better, and this observation applies to the current one as well.
Peer review completion
Journal article(s) based on this preprint
Viewed
HTML | XML | Total | BibTeX | EndNote | |
---|---|---|---|---|---|
203 | 83 | 27 | 313 | 17 | 11 |
- HTML: 203
- PDF: 83
- XML: 27
- Total: 313
- BibTeX: 17
- EndNote: 11
Viewed (geographical distribution)
Country | # | Views | % |
---|
Total: | 0 |
HTML: | 0 |
PDF: | 0 |
XML: | 0 |
- 1
Cited
1 citations as recorded by crossref.
Juliane Dannberg
Wolfgang Bangerth
Elbridge Gerry Puckett
Cedric Thieulot
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(5240 KB) - Metadata XML