Islet: Interpolation semiLagrangian elementbased transport
 Sandia National Laboratories, Albuquerque, New Mexico
 Sandia National Laboratories, Albuquerque, New Mexico
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). SemiLagrangian (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 semiLagrangian (ISL) methods contains potentially extremely efficient SL methods. We describe a set of ISL bases for elementbased 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 GaussLegendreLobatto (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 finitevolume 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 gmd2021296', 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.geoscientificmodeldevelopment.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 longterm 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

AC1: 'Reply on CEC1', Andrew Bradley, 25 Oct 2021
The comment was uploaded in the form of a supplement: https://gmd.copernicus.org/preprints/gmd2021296/gmd2021296AC1supplement.pdf

AC1: 'Reply on CEC1', Andrew Bradley, 25 Oct 2021

RC1: 'Comment on gmd2021296', Anonymous Referee #1, 29 Nov 2021
Review of "Islet: Interpolation semiLagrangian elementbased transport" by Andrew Bradley, Peter Bosler, and Oksana Guba.
The paper presents a new set of basis functions for semiLagrangian 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 semiLagrangian advection schemes, and 2) its implementation in an atmospheric model using three different (subelement) 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 3grid 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 spacetime 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 ith 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 spacetime matrix corresponds to a target node."
You should define the spacetime matrix for this statement to be comprehensible.l177: Have you defined the basis to be symmetric somewhere?
l195: "L provides a basis for degreed 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 nvector 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 ddegree 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 hhalo 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 adhoc; can you guarantee that density values are realistic?
scaling figures 68, 1116: 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 ...

AC2: 'Reply on RC1', Andrew Bradley, 07 Dec 2021
The comment was uploaded in the form of a supplement: https://gmd.copernicus.org/preprints/gmd2021296/gmd2021296AC2supplement.pdf

AC2: 'Reply on RC1', Andrew Bradley, 07 Dec 2021

RC2: 'Comment on gmd2021296', Anonymous Referee #2, 04 Dec 2021
This manuscript describes an interpolationbased semiLagrangian (SL) method for the transport problem on spectralelement (SE) domains. The SL transport schemes are widely used for multitracer 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 “massfixing” usually employed for global conservation — which may have an adverse effect for climatescale (long term) integration due to the local mass drifting. On the other hand, a finitevolume 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 (rezoning), a limiter or shapepreserving 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 highorder spectralelements are very challenging. Authors have proposed an interpolationbased SL method Islet for the SE discretization. Instead of using the “unstable” native highorder 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 massfixing. The authors argue that the Islet method can handle tracer transport as well as the remapping between physics & dynamics grids, and incorporate shapepreservation 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 twopart 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 twopart 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 115120) 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 highorder 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 SLSE 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 elementbased 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 470475). 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 multitracer 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 massfixing is the way to go, one could use the RBFbased (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 highorder elementbased 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.

AC3: 'Reply on RC2', Andrew Bradley, 07 Dec 2021
The comment was uploaded in the form of a supplement: https://gmd.copernicus.org/preprints/gmd2021296/gmd2021296AC3supplement.pdf

AC3: 'Reply on RC2', Andrew Bradley, 07 Dec 2021

EC1: 'Comment on gmd2021296', 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.

AC4: 'Reply on EC1', Andrew Bradley, 07 Dec 2021
The comment was uploaded in the form of a supplement: https://gmd.copernicus.org/preprints/gmd2021296/gmd2021296AC4supplement.pdf

AC4: 'Reply on EC1', Andrew Bradley, 07 Dec 2021

AC5: 'Comment on gmd2021296', Andrew Bradley, 15 Dec 2021
The comment was uploaded in the form of a supplement: https://gmd.copernicus.org/preprints/gmd2021296/gmd2021296AC5supplement.pdf
Andrew M. Bradley et al.
Data sets
SL GPU performance data on Summit Mar 2021 Andrew M. Bradley https://github.com/E3SMProject/perfdata/tree/main/nhxxslsummitmar2021
Model code and software
Islet methods Andrew M. Bradley, Peter A. Bosler, Oksana Guba https://github.com/E3SMProject/COMPOSE/releases/tag/v1.1.1
Performanceportable implementation in HOMME dycore Andrew M. Bradley, Peter A. Bosler, Oksana Guba https://github.com/ambrad/E3SM/releases/tag/islet2dpapersummitslgputimings
Andrew M. Bradley et al.
Viewed
HTML  XML  Total  BibTeX  EndNote  

849  195  29  1,073  12  8 
 HTML: 849
 PDF: 195
 XML: 29
 Total: 1,073
 BibTeX: 12
 EndNote: 8
Viewed (geographical distribution)
Country  #  Views  % 

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