Islet: Interpolation semi-Lagrangian element-based transport

Bradley, Andrew M.; Bosler, Peter A.; Guba, Oksana

doi:https://doi.org/10.5194/gmd-2021-296

Preprints

https://doi.org/10.5194/gmd-2021-296

Preprints

Submitted as: development and technical paper

27 Sep 2021

Submitted as: development and technical paper | 27 Sep 2021

Status: a revised version of this preprint is currently under review for the journal GMD.

Islet: Interpolation semi-Lagrangian element-based transport

Andrew M. Bradley, Peter A. Bosler, and Oksana Guba Andrew M. Bradley et al.

Show author details

Sandia National Laboratories, Albuquerque, New Mexico

Received: 27 Aug 2021 – Discussion started: 27 Sep 2021

Abstract. Advection of trace species (tracers), also called tracer transport, in models of the atmosphere and other physical domains is an important and potentially computationally expensive part of a model's dynamical core (dycore). Semi-Lagrangian (SL) advection methods are efficient because they permit a time step much larger than the advective stability limit for explicit Eulerian methods. Thus, to reduce the computational expense of tracer transport, dycores often use SL methods to advect passive tracers. The class of interpolation semi-Lagrangian (ISL) methods contains potentially extremely efficient SL methods. We describe a set of ISL bases for element-based transport, such as for use with atmosphere models discretized using the spectral element (SE) method. An ISL method that uses the natural polynomial interpolant on Gauss-Legendre-Lobatto (GLL) SE nodes of degree at least three is unstable on the test problem of periodic translational flow on a uniform element grid. We derive new alternative bases of up to order of accuracy nine that are stable on this test problem; we call these the Islet bases. Then we describe an atmosphere tracer transport method, the Islet method, that uses three grids that share an element grid: a dynamics grid supporting, for example, the GLL basis of degree three; a physics grid with a configurable number of finite-volume subcells per element; and a tracer grid supporting use of our Islet bases, with particular basis again configurable. This method provides extremely accurate tracer transport and excellent diagnostic values in a number of validation problems. We conclude with performance results that use up to 27,600 NVIDIA V100 GPUs on the Summit supercomputer.

Andrew M. Bradley et al.

Status: final response (author comments only)

CEC1:
'Comment on gmd-2021-296', Juan Antonio Añel, 25 Oct 2021

Dear authors,

After checking your manuscript, it has come to our attention that it does not comply with our Code and Data Policy.

https://www.geoscientific-model-development.net/policies/code_and_data_policy.html

You have archived your code in GitHub. However, GitHub is not a suitable repository. GitHub itself instructs authors to use other alternatives for long-term archival and publishing, such as Zenodo. Therefore, please, publish your code in one of the appropriate repositories. These repositories must be available for the review process, as preliminary versions are equally relevant as final ones.

You must reply to this comment, including the link to the repositories and DOIs.

Also, you must include in a potential reviewed version of your manuscript the modified 'Code and Data Availability' section.

Please, note that in the current version of the manuscript, one of the hyperlinks to one of the GitHub repositories is broken.
Juan A. Añel
Geosc. Mod. Dev. Executive Editor

Citation: https://doi.org/10.5194/gmd-2021-296-CEC1
- AC1: 'Reply on CEC1', Andrew Bradley, 25 Oct 2021
  
  The comment was uploaded in the form of a supplement: https://gmd.copernicus.org/preprints/gmd-2021-296/gmd-2021-296-AC1-supplement.pdf
  
  Citation: https://doi.org/10.5194/gmd-2021-296-AC1
RC1:
'Comment on gmd-2021-296', Anonymous Referee #1, 29 Nov 2021

Review of "Islet: Interpolation semi-Lagrangian element-based transport" by Andrew Bradley, Peter Bosler, and Oksana Guba.

The paper presents a new set of basis functions for semi-Lagrangian advection method in spectral element models. The presented method is novel and its performance in atmospheric test cases is encouraging. The paper therefore warrants to be published. The main drawback of the paper is that it is quite tedious to read. This paper essentially presents two things: 1) a novel optmized basis functions for 1D semi-Lagrangian advection schemes, and 2) its implementation in an atmospheric model using three different (sub-element) grids. Both of these are quite complex topics, and their discussion is intertwined (e.g. in section 2) and the reader is easily lost in details. The overall presentation should be improved before the paper can be accepted for publication. Considering the amount of work, I recommend a major review.

Major comments:

For clarity, I suggest that you clearly define the two considered problems, 1D SL advection method, and the implementation in 3-grid atmospheric models, already in the introduction. The introduction is now quite short and actually does not mention many important concepts relevant for the paper.

Section 2 is rather difficult to follow, consider revising. I suggest to start by defining the 1D discretization with N elements, the interpolant functions within each element, and the properties of the interpolant functions (e.g. basis functions, continuity and symmetry). Presently, the interpolant functions first appear only in section 3.1. A figure could clarify the concepts, including the source and target elements/nodes. I would also define a symbol for the basis functions themselves, instead of using L from section 3.1 (L is the interpolant function itself). The discussion of stability becomes comprehensible only after the discretization has been introduced.

As I understand, the 3 axioms of advection methods are: 1) global conservation, 2) preservation of constant tracers (sometimes called local conservation, or tracer consistency), and 3) monotonicity (i.e. no spurious overshoots appear). These properties apply to both the SL tracer advection scheme and the remap operators between grids. In section 1.2 these concepts seem to be mixed and referred to by different names (plausibly due to historical reasons; the so called "property preservation" is just a combination of 1, 2, and 3). Consider revising.

Figures 20 and 21 shows larger diffusion for shorter time steps. Could you elaborate on this? Is the traditional SL property that longer time steps reduce diffusion (while the solution can degrade in other metrics)? How would one choose the right time step in practice? In the case of Eulerian transport it is easy: take the maximum stable one, and you are guaranteed to satisfy all the necessary properties.

Minor comments:

Throughout the manuscript the authors use the terms "mixing ratio" and "tracer" for the advected quantity q_i, seemingly interchangeably. For the sake of clarity I would prefer just to use "tracer".

l35: what is the definition of "local conservation" here?

l155: "This structure arises as follows. Consider a continuous discretization using a nodal n_p -basis, n_p = d+1, with n_p the number

of nodes. The grid has N elements. Each row of the space-time matrix corresponds to a target node."

This description is too brief to be understandable, please elaborate. This is my interpretation of the discretization: The 1D domain is divided into N elements. The solution in each element is approximated by a continuous function, defined by n_p basis functions. Thus a function f in an element e can be written as f_e(x) = sum_i=1^{n_p} f_i psi_i(x) where psi_i and f_i denote the i-th basis function and its corresponding coefficient. Each basis function is associated with a node x_i within the element; The basis is Lagrangian (a.k.a. nodal), i.e. f_i(x_i) = 1. Furthermore, the discretization of the function is continous across element interfaces, implying that the neighboring elements share a (exactly one?) basis function. Furthermore, the basis is assumed to be symmetric about the center point of the element.

l156: "Each row of the space-time matrix corresponds to a target node."

You should define the space-time matrix for this statement to be comprehensible.

l177: Have you defined the basis to be symmetric somewhere?

l195: "L provides a basis for degree-d polynomials." I would say that the basis functions are the \Pi_i functions defined in the equation of L; L itself is the interpolant function defined by the basis and the specific nodal values y(i).

l195: "These are supported by n = d + 1 points, each an element in the n-vector xn." To be consistent with the literature I would use the term "node" instead of "point".

l204: "Given a departure point x" These properties define the interpolant functions, thus there's no need to say that x is a departure point, it can be any point within the element.

l208: I think this constraint is equivalent to saying that the basis must Lagrangian or nodal?

l246: only here you define a d-degree polynomial. This would be useful already in section 2.

l257: what is Runge's phenomenon? what is Lebesque constant? help the reader to understand the rationale behind your work.

Section 3.6: the description of the search algorithm is quite technical and could perhaps be moved to the appendix; it is not necessary to follow the main storyline of the paper.

Section 4: mention TTPR already in the introduction as it seems to be relevant for the entire Islet method.

l439: earlier tracer tendencies were denoted by f_i \Delta t

l442: "immediate element neighbors" Are these neighbors that share an edge or vertex?

l464: "In contrast, ..." Meaning unclear, please revise.

l507: What is "tracer density"? Is it just \rho? Then, for clarity, I'd call it "density"

section 4.3: while computational efficiency is important, I would move this section to the appendix, as it

l570: "is proportional to the number of grid points" Should read: "proportional to the square of number of grid points"

l572: what is a naive h-halo exchange and how does it differ from what is proposed here?

section5, l606: please mention the problem domain (full sphere?) and the equations that are being solved (pure advection on the tracer grid?)

l631: why can you not use the dynamical core to compute the density? the chosen approach seems rather ad-hoc; can you guarantee that density values are realistic?

scaling figures 6-8, 11-16: for easier readability do not use red line twice for np=6 and np=12. x axis label is missing.

l703: "has even more accuracy" Can you quantify this? Is the difference significant?

l711: "For each value ..." unclear sentence, please revise.

l741: 0 in m(0) stands for time t=0?

l780: what is a terminator?

l897: Here you define the basis functions \phi_i. This notation would be useful throughout the manuscript, already in Sect 2 and 3.1.

abstract and intro: you define an abbreviation "dycore". I would omit it as it does not really save space, and it is not used frequently in the paper.

Typos:

l745: "shows the results"

l815: SYPD numbers printed above ...

Citation: https://doi.org/10.5194/gmd-2021-296-RC1
- AC2: 'Reply on RC1', Andrew Bradley, 07 Dec 2021
  
  The comment was uploaded in the form of a supplement: https://gmd.copernicus.org/preprints/gmd-2021-296/gmd-2021-296-AC2-supplement.pdf
  
  Citation: https://doi.org/10.5194/gmd-2021-296-AC2
RC2:
'Comment on gmd-2021-296', Anonymous Referee #2, 04 Dec 2021

This manuscript  describes an interpolation-based semi-Lagrangian (SL) method for the transport problem on spectral-element (SE) domains.  The SL  transport  schemes are  widely used  for multi-tracer transport  in atmospheric models due to their  accuracy and computational efficiency.  The classical SL method employs interpolation at the upstream locations of the backward trajectories to estimate the advecting scalar values at the new time level. However, such an approach is not conservative per se, for practical applications an arbitrary procedure  known as the “mass-fixing” usually employed for global conservation — which may have an adverse effect for climate-scale (long term) integration due to the local mass drifting. On the other hand, a finite-volume formulation of  the SL method is conservative by design, where the upstream interpolation over the Lagrangian element is replaced by integration constrained to be locally  (hence globally) mass conservative. The conservative data transfer from regular Eulerian grid to the deformed Lagrangian grid  often referred to as the remapping (re-zoning),  a limiter or shape-preserving scheme is usually employed for physically realizable solutions.  A wide body of literature is available for both conservative and classical SL methods.

Implementation of  conservative SL method on spherical domains tiled with high-order spectral-elements are very challenging. Authors have proposed an interpolation-based SL method  for the SE discretization.  Instead of using the “unstable” native high-order  interpolator (basis function) they have devised a cumbersome numerical procedure   which employs an alternative grid system within each  spectral element, adding another layer of complexity. The Islet method is not conservative, nevertheless, the global conservation is achieved by mass-fixing. The authors argue that the Islet method can handle tracer transport as well as  the remapping between   grids, and incorporate shape-preservation filters.

The manuscript is very long, the Islet interpolation as described by the authors is extremely complex. Authors failed to explain the core interpolation algorithm with clarity, there are many statements in the manuscript which leads to ambiguity.  The numerical analysis part is very intense  maybe more suitable for a  computational math journal  (e.g., SIAM / JCP)  than  the  GMD. The subject covered could be split into a two-part paper, one describing the basic algorithm and analysis with more details and rigor, and the second part for implementation and validation with standard tests. This would be helpful for better reading. The current  manuscript is written in an awkward manner  and is unacceptable for publication.

Recommendation:      Major revision, possibly resubmit  as a two-part manuscript. Authors should address the following questions.

(1) The stability associated with the SL method is that the deformational Courant number  (Lipschitz condition) should not exceed unity, in plain  language, the trajectories should not cross intersect (see, Staniforth & Cote’s 1992 MWR paper). Is the  cubic ISL method (lines  115-120) unstable due to this condition? Need some explanation.

(2) The SL  transport scheme can be stabilized using a limiter,  filter  or with an explicit diffusion (see, Ullrich & Norman, QJRMS, 2014). You can use the native high-order SE interpolation (basis function) for the SL transport combined with the limiter which you are already using for the Islet method. It will be interesting to see how the Islet method compares with this simple SL-SE scheme employing 4x4 GLL grid (I guess that is the SE grid choice made for the operational E3SM).

(3) It is not convincing to have 3 grid systems (physics: FV, dynamics: GLL, transport: tweaked GLL)  in a SE modeling framework. The Fig.5  shows such a grid configuration, and it appears to be very  challenging. At  a very high (NH)  resolution the data movement is a major issue for an element-based  Galerkin model (DG/SE).   A typical climate model may have O(100) tracers, an additional   tracer grid with more  DOF than the dynamic grid can exacerbate this problem. This will limit the use of Islet scheme, how do you address it?

(4) With real data you have velocity information only available at the  GLL (dynamics) grid, the way you find   the 2D trajectory information   using the 3D Cartesian coordinates  leads to additional computational overhead when the method is extended to the 3D application (line 470-475). This needs some justification, why not use the spherical (u,v) components or corresponding contravariant vectors?  It is not clear that the maximum eigenvalue required for the interpolation is the tracer data dependent, in that case you have a serious  computational overhead for the  multi-tracer applications, Please clarify!   What is the computational halo  requirement for an  SE  stencil with  NxN GLL points, when the shape preserving limiter is applied?

(5) What is the special advantage of using  Islet  method? It seems you have introduced a complex numerical method for a relatively simple linear transport problem.  If mass-fixing is the way to go, one could use the RBF-based (Kriging type)  interpolator which provides very accurate solution, and no need for  the expensive search for max eigenvalue etc.

(6)  The results are looking good, authors should limit the number of figures and make an effort to compare the results  with  that of other high-order element-based schemes. Why the results from your own previous papers  (Bosler et. al. 2019, SIAM J Sci. Computing; Guba et al. 2014, JCP) discussed? These results  should be compared and  the relative merits should be discussed.

(7) There are many undefined terms (e.g. CAAS) and notations which I am going to list, this should be fixed.

Citation: https://doi.org/10.5194/gmd-2021-296-RC2
- AC3: 'Reply on RC2', Andrew Bradley, 07 Dec 2021
  
  The comment was uploaded in the form of a supplement: https://gmd.copernicus.org/preprints/gmd-2021-296/gmd-2021-296-AC3-supplement.pdf
  
  Citation: https://doi.org/10.5194/gmd-2021-296-AC3
EC1:
'Comment on gmd-2021-296', James Kelly, 04 Dec 2021

In light of the second reviewer's comments, I am encouraging the authors to withdraw the manuscript and submit to another journal, such of Journal of Comptuational Physics.

Citation: https://doi.org/10.5194/gmd-2021-296-EC1
- AC4: 'Reply on EC1', Andrew Bradley, 07 Dec 2021
  
  The comment was uploaded in the form of a supplement: https://gmd.copernicus.org/preprints/gmd-2021-296/gmd-2021-296-AC4-supplement.pdf
  
  Citation: https://doi.org/10.5194/gmd-2021-296-AC4
AC5: 'Comment on gmd-2021-296', Andrew Bradley, 15 Dec 2021

The comment was uploaded in the form of a supplement: https://gmd.copernicus.org/preprints/gmd-2021-296/gmd-2021-296-AC5-supplement.pdf

Citation: https://doi.org/10.5194/gmd-2021-296-AC5

Andrew M. Bradley et al.

Data sets

SL GPU performance data on Summit Mar 2021 Andrew M. Bradley https://github.com/E3SM-Project/perf-data/tree/main/nhxx-sl-summit-mar2021

Model code and software

Islet methods Andrew M. Bradley, Peter A. Bosler, Oksana Guba https://github.com/E3SM-Project/COMPOSE/releases/tag/v1.1.1

Performance-portable implementation in HOMME dycore Andrew M. Bradley, Peter A. Bosler, Oksana Guba https://github.com/ambrad/E3SM/releases/tag/islet-2d-paper-summit-sl-gpu-timings

Andrew M. Bradley et al.

Viewed

Total article views: 1,073 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
849	195	29	1,073	12	8

HTML: 849
PDF: 195
XML: 29
Total: 1,073
BibTeX: 12
EndNote: 8

Views and downloads (calculated since 27 Sep 2021)

Month	HTML	PDF	XML	Total
Sep 2021	89	22	4	115
Oct 2021	145	31	4	180
Nov 2021	107	19	1	127
Dec 2021	160	46	14	220
Jan 2022	86	18	1	105
Feb 2022	80	12	2	94
Mar 2022	74	15	2	91
Apr 2022	68	24	0	92
May 2022	40	8	1	49

Cumulative views and downloads (calculated since 27 Sep 2021)

Month	HTML	PDF	XML	Total
Sep 2021	89	22	4	115
Oct 2021	145	31	4	180
Nov 2021	107	19	1	127
Dec 2021	160	46	14	220
Jan 2022	86	18	1	105
Feb 2022	80	12	2	94
Mar 2022	74	15	2	91
Apr 2022	68	24	0	92
May 2022	40	8	1	49

Viewed (geographical distribution)

Total article views: 1,014 (including HTML, PDF, and XML) Thereof 1,014 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 20 May 2022

Short summary

Tracer transport in atmosphere models can be computationally expensive. We describe a flexible and efficient interpolation semi-Lagrangian method, Islet. It permits using up to three grids that share an element grid: a dynamics grid for computing quantities such as the wind velocity; a physics parameterizations grid; and a tracer grid. Islet performs well on a number of validation problems and achieves high performance on up to 27,600 NVIDIA V100 GPUs on the Summit supercomputer.


Total:	0
HTML:	0
PDF:	0
XML:	0