the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 08 Oct 2024
Anisotropic metric-based mesh adaptation for ice flow modelling in Firedrake
Abstract. Glaciological modelling is a computationally challenging task due to its high cost and complexity associated with large spatial- and long time-scale simulations. In this paper, we provide a comprehensive overview of state-of-the-art feature-based anisotropic mesh adaptation methods and demonstrate their effectiveness for time-dependent glaciological modelling using the Python-based Firedrake finite element library. We introduce a novel hybrid time-dependent fixed-point mesh adaptation algorithm that generates a more optimal initial mesh sequence. The algorithm requires approximately 50 % fewer iterations in order to reach mesh convergence, while still controlling spatial error and its temporal distribution. We demonstrate the effectiveness of anisotropic mesh adaptation and the novel fixed-point algorithm on a Marine Ice Sheet Model Intercomparison Project (MISMIP+) experiment. We show that we are able to achieve solution accuracy comparable to a uniform 0.5 km resolution mesh simulations by using a sequence of adapted meshes with, on average, 10–30 times fewer vertices, depending on the sensor field used to drive mesh adaptation. Due primarily to the iterative nature of the mesh adaptation process employed, this translates in practice into a 3–6 times lower overall computational cost.
- Preprint
(2668 KB) - Metadata XML
- BibTeX
- EndNote
Status: open (until 13 Jan 2025)
-
RC1: 'Comment on egusphere-2024-2649', Ed Bueler, 05 Dec 2024
reply
Summary: This manuscript describes, implements, and tests an adaptive mesh refinement (AMR) scheme for the 2D time-dependent shallow shelf approximation (SSA) models of marine ice sheets. The scheme is based on metric-based mesh adaptation, using fields derived a posteriori from ice thickness, ice flow speed, the magnitude of basal shear stress, or intersections thereof, as the local refinement criterion. The essential idea is that the adapted mesh comes from solving an optimization problem for a Riemannian metric field, constrained by a fixed number of mesh nodes, using Hessians (2nd derivatives) of the current numerical solution. The scheme has an open-source Firedrake-based Python implementation, with mesh construction from the computed metric ("mesh adaptation") via the library Mmg2d. The AMR scheme is tested on the MISMIP+ marine ice sheet idealized experiment. With the caveat that no MISMIP+ exact solution is known, clear reduction in mesh complexity is observed, relative to uniform mesh refinement, for a given numerical error level. Modest speedup, in a time-to-error sense, occurs in the right circumstances, but the large cost of computing solution Hessians is noted.
Recommendation: The ideas here, and I believe the implementation (though I have not tested it), represent valuable work. Concretely, an open source, Firedrake-based toolchain for adding principled (metric and anisotropic) AMR to marine ice sheet simulations is a *good thing*. The reported performance results suggest a meaningful speed-up over a uniform mesh, though apparently less than an order of magnitude. (In this regard, speedup must at least compensate for the increased complexity.) This work should be published in a well-written and self-contained form. However, while the authors intend that the current manuscript is "self-contained", in fact it is fundamentally lacking in this regard. Almost all readers will have difficulty in following the technical details; only those with prior metric-based adaptation experience could be comfortable. Serious readers will split into those who willing use the authors' implementation, trusting that things work in some reasonable manner, and a (smaller?) group who must reconstruct for themselves, from the literature and the source code, many aspects of how things actually work, either in high-level outline or in detail. I care about the latter group; please see the suggestions following. I would suggest that the paper be returned for major revisions.
Suggested presentation, so as to be self-contained: The content in sections 2.1--2.4 is reasonably self-contained, but it would be more so if it explicitly described how things work for the familiar Poisson equation, presumably with the solution (absolute value of the solution?) as the "sensor field". A picture of how a mesh is actually generated from a metric is currently missing as far as I can tell, and it could be easily added to this Poisson case; it is the main point. The fact that actual mesh adaptation occurs inside Mmg2d should be made clearer to the casual reader. A clear visual image of metric and mesh anisotropy is also missing, and this could be added. Now a new section 3 could be created containing the marine ice sheet equations (which should be stated as a coherent, essentially complete block), the new time-dependent ideas, and a much clearer discussion of the nontrivial sensor field choice for marine ice sheets. These new things are the novel ones, relative to prior anisotropic metric-based AMR, as far as I can tell. Such important new concepts should not be confusingly scattered in the way they currently are. It must be clearly indicated to the glaciologically-inclined reader, presumably in the new section 3, and in any case well *before* the end of the results section, that vector-valued SSA-based marine ice sheet simulation requires a nontrivial physics choice regarding AMR, and that this choice is not a settled matter. That is, one needs to *choose* between surface speed, ice thickness, basal stress, intersected speed/thickness, and etc as the sensor field. (The current manuscript assumes glaciological readers will happily read all the way to page 20 and 3.4.2 before having any clear statement of the physics basis on which the mesh is refined! This is crazy.) The same choice should be clearly sketched in the Introduction, and perhaps even the abstract. Then, in section 4 on results, the particular MISMIP+ experimental configuration could be described, and finally the actual results given. Finally, many of the Figures are elaborate in the modern Nature/Science style, and thus no reader understands them; Figures 5,9,10 are this way and I suggest they be broken into discrete digestable form. The spaghetti needs to be unwound.
Specific Comments on Manuscript
line 2,3 (in abstract): Avoiding the boring sales pitch about "comprehensive overview of state-of-the-art" stuff. Perhaps "In this paper we propose feature-based anisotropic mesh adaptation ..."
line 6: Remove "more optimal"; it does not mean anything. In the same line, "50% fewer" than what? Say briefly.
line 11: Almost no readers (me included) know what a "sensor field" is. This needs a rewrite in the abstract. I suggest avoiding "sensor field" here (but paraphrase it), and then carefully define it later, e.g. at line 181, which seems to be the next use of the phrase.
line 11: What does "Due primarily to the iterative nature" mean here? The reader is unlikely to believe that just because it is iterative it is therefore faster or better. Perhaps this last sentence should be something like: "In practice the adapted mesh technique translates into a ..."
lines 14-81: I think the Introduction is quite well-written, and it accords with my understanding of the development of models suitable for marine ice sheets.
line 22: For the reader to understand the meaning here I suggest: "... the majority have focussed on reducing computational costs by numerical and meshing techniques, while maintaining the model physics."
line 27: Recent significant progress has been made on this contact problem. See: de Diego, G. G., Farrell, P. E., & Hewitt, I. J. (2022). Numerical approximation of viscous contact problems applied to glacial sliding. Journal of Fluid Mechanics, 938, A21.
line 29: It is not clear what "Alternative" means here. This lead has the same meaning?: "Mesh adaptation methods select finer ..."
line 65: Suggest "allows to equidistribute" --> "equidistributes".
line 72: "velocity norm" --> "speed"
lines 73--75: I had to reread the last two sentence of this paragraph several times. I suggest this 3 sentence replacement: "These methods generate only a single mesh that is used throughout the simulation. However, they do not consider the nonlinear and time-dependent coupling between the solution and the underlying mesh. This coupling suggests an iterative mesh adaptation procedure, which we will implement."
line 76: As noted earlier, my main concern about the current manuscript is that it is *not* a sufficiently self-contained description, given the reader audience, including myself.
line 79: "but implement" --> "and we implement"
line 83: Presumably \Omega has a polygonal boundary.
line 88-89: It is true that there is not a "monotonic" relationship, but that is hardly the point. AMR does not guarantee monotonicity here either, as far as I know. The point is that uniform refinement can lead to only small, thus expensive, reductions in the error.
equation (1), page 4: The displayed equation should be taken more seriously. It should say something like "find \H_opt solving min E(\H) over \H with N_v vertices". Then the sentence after can be shortened.
lines 96-97: Suggest: "In anisotropic mesh adaptation the optimal mesh is only found approximately, and by iteration."
line 102: "a way" --> "a local way"
lines 103-104: Actually put the definition after the sentence which promises the definition: "A Riemannian metric is a smooth function yielding a positive-definite 2x2 matrix M(x) at each x in the domain."
somewhere around equation (2): State that \lambda_1 \ge \lambda_2 somewhere, as I am pretty sure you are using that.
line 115: "a metric tensor M is a continuous analogue of a single mesh element" makes no literal sense. Probably: "a metric tensor field M(x) is the continous analog of the changing shape [dimensions?] of elements as one moves through the mesh." Yes? Otherwise "continuous" is not modifying anything?
lines 119-120: "fully" makes no sense here except as vacuous advertising. How about "... in R), thus the metric-based approach allows anisotropy in the resulting mesh."
caption Figure 1: "is 16 times" --> "is in fact 16 times"
caption Figure 1: ", which is not reflected in the figure for readability" --> "; the domain is laterally shortened for legibility."
line 121: "metric space" actually means something, which is not here. How about "... the Riemannian space \M = (\Omega, M(x)) models ..." Surely one can call M(x) a metric tensor field, or just a metric field, and not conflate the ideas of "metric space" and "Riemannian manifold"?
line 126: "is of" --> "has"?
line 127: Since I actually know what a Riemannian manifold is, and what a metric space is, I am again irritated that you could say "Riemanninan space" or "Riemannian manifold" here without harm, but instead you have a jumble. The metric space structure of a Riemannian space/manifold comes from minimizing the integral of the Riemannian metric over curves, right? But you don't need that metric space structure, I am quite sure. In conclusion, avoid "metric space" unless that is really what you mean.
following up: I see that the fundamental confusion is already present in Alauzet 2010. That author simply decides to measure distances between points in a "Riemannian metric space" in a way which is different from the metric space structure of a Riemannian space. He simply decides that you integrate over Euclidean lines. Indeed he does not need the correct metric space, and just overwrites it with his stuff using the same words. It is sad and destructive that that author has now generated a branch of mathematics where words mean different things from the rest of math. (You will never be able to do metric-based adaptation on a Riemannian manifold without undoing his language choice.) Why did he call his thing "Riemannian" when he announces that it is not? I would be cool with "Alauzet metric space".
line 135: Just strike "which appear to be ...". The shortening should be clearly stated in the Figure caption, e.g. as suggested above. Then do not refight that battle; it is confusing to do so.
line 148: ||u - u_h|| is *not* the "local interpolation error". Correct this; I think you are merely talking about ||u-\Pi_h u||?
end of line 153: "interpolant", right? (I see Alauzet already makes this English mistake.)
equation (6): How about
"find \M_{L^p} which minimizes \E_{L^p}(\M) = \int_\Omega ..."
That is, the input variable to the objective is (caligraphic) \M. Putting the L^p subscript on "\M_{L^p}" then labels this as the minimizer of \E_{L^p}.line 164: What does "smooth" mean? Drop that word? This is an L^2 projection into a continuous FE space, right?
around lines 155-164: I think it should be said that you don't actually have access to u and you are applying the idea to u_h a posteriori.
line 166: At this point the reader has no idea what you will later mean by "one partition of the simulation time interval". Just refer to a "time interval" here, and don't burden the reader with your later optimizations.
line 176: "may necessitate" --> "needs"
line 181: At some point you need to define "sensor field". As done here it is either too late or too awkwardly said. Maybe connect to equation (6) at the beginning of this paragraph; say that in the scalar Poisson context one uses u for adaptation, and one calls this use the sensor field. Now say that there are choices in the time-dependent coupled context.
line 184: "mean" --> "method"
equation (10): Here I think this is only well defined if lambda_1 >= lambda_2, right? See comment about equation (2).
line 190: "eigenvalues" --> "eigenvalue estimates"? I don't see how they can be eigenvalues unless the two metrics are assumed to be simultaneously diagonalizable by the {e_i}?
lines 191, 192: strike "popular" and "commonly"
Figure 2: Note this is a detailed figure about a fiddly detail. It should go later, as part of putting novel ideas later. A couple of early figures emphasizing the main ideas of metric-based anisotropic mesh adaptation should be substituted.
section 2.4: This immediate descent into details of postprocessing misses the main message: from the a posteriori computed metric estimate M(x), or M^i(x) in the time-dependent case, you are asking Mmg2d to generate mesh(es). *Then* talk post processing details, though probably put later.
line 238: Speaking as a reader, I would want the "several shortcomings" to be explained here. Evidently I am too stupid to remember the brilliant exposition from the Introduction.
line 245: Interleaving the explanation of the algorithm and advertising its benefits is not helpful. Do the former first.
line 247: Probably clearer without "before they are adapted".
line 248: Again, advertising is mixed into algorithm description. Not helpful.
line 252: On first pass I definitely did not understand the meaning of "The fixed-point algorithm iterates ..." *here*. (I wrote: """Is this recomputing solutions u along the way? Are you returning to the start of the time partitions to revise the metric, or is the metric getting updated only in a feed-forward manner. Iteration is an ambiguous word, but your use of "fixed-point" suggests things are not feed-forward only. Does the "x" symbol in Figure 2 mean "compute a new metric but don't change the mesh"?""")
section 2.5 generally: The first 4 paragraphs are quite good. Starting in the 5th paragraph the job is to describe unambiguously an algorithm. A pseudocode may be better than the paragraph-with-interleaved-commentary-and-advertising current form. Figures showing the action on coarse meshes might help.
generally: These time-dependent mesh adaptation schemes seem to be awkward, and early-development. Do you really recommend them for production simulations?
section 2.6: This part makes sense to me, but how understandable is it if you are not already comfortable with supermesh ideas, interpolation, and Galerkin projection?
line 290: The avoidance of negative values via P^1 makes sense, but I suppose this begs the question of whether one should, instead of h-refinement and the costs of AMR, instead p-refine in a way which avoids the out-of-bounds problem. See recent work: Kirby & Shapero 2023. High-order bounds-satisfying approximation of partial differential equations via finite element variational inequalities arXiv:2311.05880
line 305: Give some citation for "PETSc's Riemannian metric functionality"?
line 308: As a nonexpert I don't know what "mesh optimization" means here. Which of the many ideas relating to building better meshes does it refer to?
sections 2.6, 2.7: Easy to read.
line 340: "the equations of glacier flow" should actually be stated. (The "the" is ridiculous here; it is a field with lots of models.)
line 342: Indeed this alternation is "as with most current ... models". The scheme is called the explicit Euler approximation to the coupled model, that is, the scheme is explicit and first-order! As with current ice sheet models, attaching the "prognostic" and "diagnostic" language, a leftover of pre-numeric weather forecasting, is not helpful. The current authors are not to blame here, but they don't have to climb into the dumpster too.
line 353: Arg. Is "P" the pressure or the effective pressure? Effective pressure is the *difference* between the ice pressure (overburden pressure) and the pressure of subglacial water. At the bottom of floating ice this value is zero; assuming h is ice thickness and d_W is the depth below the water line of the bottom of the ice, then rho_i g h = rho_W g d_W because of flotation. Please don't use "P" for the effective pressure. Fix all this.
line 356: Please do not assert that the scheme in IcePack is unconditionally stable. It is not. Importantly, Shapero (2021) does not assert that it is. They use a scheme which is unconditionally stable *for the linear advection equation*. Then they say: "We note that, while the stability properties of different schemes for the linear advection equation guided our choice of method, the coupled system for both thickness and velocity is not linear and not hyperbolic." The scheme in IcePack is implicit for the linear advection equation but *not* implicit for the coupled system. In this context it completely suffices to say that "Finally, IcePack implements a second-order Lax-Wendroff scheme for transport (Shapero, 2021)."
Figure 4: If this result is worth showing for uniform refinement, then also show it for a preferred AMR strategy? Presumably, with fewer dofs than the 0.25km uniform?
line 363-364: Me too. You don't need to include this detail particularly, but the cp linesearch in PETSc never seems to work for me, and I fall back to backtracking.
line 390: "emulate a more general problem". Huh? This sentence probably means "In the mesh adaptation experiments to follow, neither exact nor high-fidelity solutions are assumed to be available." On the other hand, a couple of sentences later you treat the 0.25km uniform result as high-fidelity. Rewrite this paragraph?
line 396: I don't know what "effectively ... unbounded" means, just after you give bounds. Frankly, I have a sense that something is numerically unstable etc., and you are special pleading.
section 3.3 generally: I think you can write up your experiment choices more declaratively and less defensively. To first approximation, just say what you did.
Figure 5: A busy figure trying to show several different things, and with different horizontal axes. (Likewise Fig 9, Fig 10.) Decide to show one thing and show it?
section 3.4: Use "H" for ice thickness. "h" already has its usual meshing meaning of mesh size.
line 431: This casual mention of using \tau_b for mesh adaptation seems to be the first contact with glaciological interests. On page 17.
subsection 3.4.1: Speaking as a reader who is making an effort, I am a bit lost at this point on what the difference between global and hybrid is. I see hybrid is better (= closer to uniform ref). Is hybrid more expensive? Generally much of this subsection reads as chain of thoughts of the numerical experimenter, not a purposeful statement of results.
line 442: "equal observations" means what? Re-word.
subsection 3.4.2: I have already made the point that this physics choice, of which field to adapt upon, should be put much earlier. Needless to say I am not trying to make this point again.
equations (13), (14): Also an equation for \tau_b, right? (If not I am definitely missing something ... and it is your fault.)
line 500-501: Since basal stress nulls as you cross the grounding line, perhaps intersecting tau_b and the strain rate magnitude as sensor fields? The latter is telling you about the parts of ice shelves which are interesting. Maybe the Hessian of speed is all you need, not strain rate magnitude (i.e. as sensor field); so intersect tau_b and |u|? The thickness will follow? In any case, the equations of a glacier are (at core) balancing stresses and strain rates to produce smooth velocity and thickness fields.
Figure 8: This dark dense figure tells me little and makes me want to look elsewhere. Maybe just reverse the shading.
line 503: Again, "not reflected in Fig 7 to not hinder readability" ... is hindering readability. Since this can be clarified already for Figure 1, I don't think you need to say anything here.
Figure 9 caption: State here that errors e_h, e_u are relative to 0.25km uniform ref. The reader jumping to this Figure needs that knowledge first.
lines 534-537: I appreciate this honest assessment of the technique.
line 561: "utilized libraries" --> "libraries used here"; redundant statement of Firedrake usage
lines 566-568: Mention this fact, that basal shear stress is the preferred sensor field for efficiency (error vs time) in abstract? It is actually a result, not just a "we did some simulations" statement.
line 570: "permille" is not widely used in English literature? But perfectly clear. Interesting.
lines 573-574: suggest replacing "primarily due to the iterature procedure of mesh adaptation algorithms" with "reflecting the cost of our iterative mesh adaptation procedure"
line 575: "utilising" --> "using"
line 576: "utilised" --> "used"; suggest replacing "features of metric-based adaptation" with "available metric-base adaptation techniques"
line 582: "utilised" --> "used" ... and etc ... you get the point
line 589: "necessitated" --> "required"
lines 597-599: My conclusion has been the opposite, and it is driven by experience switching from idealized experiements with smooth bed topography to modeling whole ice sheets in the BedMachine era. In the later case one has bumpy bed essentially everywhere, and bed topography always locally influences basal shear stress through thermodynamic and liquid water means. Thus everywhere needs about comparable resolution. For example, as we get good bed elevation observations in East Antarctica, then (at least in a 2D sense) everywhere will need good res., with only a slightly enhanced need for resolution near e.g. grounding lines and shear margins. Said a different way, mesh adaptation is effective when part of the domain is easy, and I don't think that is meaningfully true, except for the small area of ice shelves, for real marine ice sheets with real beds. The contrast with mesh adaptation for fluid problems within engineered geometries is strong.
line 600: My favorite line of the paper.
Citation: https://doi.org/10.5194/egusphere-2024-2649-RC1
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 179 | 62 | 33 | 274 | 4 | 6 |
- HTML: 179
- PDF: 62
- XML: 33
- Total: 274
- BibTeX: 4
- EndNote: 6
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1