MQGeometry-1.0: a multi-layer quasi-geostrophic solver on non-rectangular geometries

Thiry, Louis; Li, Long; Roullet, Guillaume; Mémin, Etienne

doi:https://doi.org/10.5194/gmd-17-1749-2024

Articles | Volume 17, issue 4

https://doi.org/10.5194/gmd-17-1749-2024

© Author(s) 2024. This work is distributed under
the Creative Commons Attribution 4.0 License.

https://doi.org/10.5194/gmd-17-1749-2024

© Author(s) 2024. This work is distributed under
the Creative Commons Attribution 4.0 License.

Articles | Volume 17, issue 4

Model description paper

|

28 Feb 2024

Model description paper |

| 28 Feb 2024

MQGeometry-1.0: a multi-layer quasi-geostrophic solver on non-rectangular geometries

Louis Thiry, Long Li, Guillaume Roullet, and Etienne Mémin

Download

Final revised paper (published on 28 Feb 2024)
Preprint (discussion started on 27 Sep 2023)

Interactive discussion

Status: closed

RC1:
'Comment on egusphere-2023-1715', Anonymous Referee #1, 06 Oct 2023

The paper presents a new implementation of a numerical solver for the multi-layer QG equations in non-rectangular domains. The numerical method is well presented and includes all relevant citations. I believe the proposed implemenation is of potential interest for the community. Overall, the paper is well-written with only a few typographical errors. I have nonetheless a few suggestions for possible improvements. The main one is that similarly to what the authors do for the barotropically unstable shielded vortex, the authors show the evolution of the enstrophy for the second numerical experiment (the near inviscid vortex in a box with a internal wall). Other minor suggestions are listed below
1) l.15 : typo "salinity"
2) bottom of page 1: add punctuation to your equations (a comma after the first one and the full stop after the second). Check all equations
3) the implemenation does not seem to include bottom topography. Can the authors make a comment on possible implementation?
4) line 36: insert `the' before `QG equations'
5) line 51 `This ad-hoc...': a citation is required.
6) line 54: I'd suggest to use `choice' over `decision' (check throughout)
7) line 57: I'd suggest to replace warrants with guarantees (which seems more natural in the context)
8) line 71 typo `machine learning'
9) line 80 remove 'us' after `leads'
10) line 85 `lives' -> 'lies' or `is computed'
11) line 90: `friction'? : The equations given are inviscid.
12) lines 93-95 I cannot understand the sentence `Moreover...'
13) line 114 `smoothness indicator' is too vague. The reader needs a little more explanation
14) l 118-119: Please keep the same word order: TVD Runge-Kutta or Rungue-Kutta TVD
15) What is the point of the last sentence at the bottom of page 10?
16) line 240: `With' -> `with'
17) line 261: `The experiment is integrated...' -> `The equations are integrated...' (an experiment is not integrated, the equations are)
18) line 292: f-plane -> $f$-plane
19) line 302: remove the unnecessary `very'
20) line 310: grammar: `configurations are' or `configuration is'
21) line 327: `on Bottom' -> `on the bottom'

Citation: https://doi.org/10.5194/egusphere-2023-1715-RC1
- AC1: 'Reply on RC1', Louis Thiry, 09 Nov 2023
  
  We would like to thank Referee #1 for their comments, suggestions, and precise list of identified typos that we will address.
  
  As recommended by the referee, we will include in the revised version a plot illustrating the evolution of the enstrophy for the vortex wall experiment.
  
  Concerning other suggestions:
  
    3) the implementation does not seem to include bottom topography. Can the authors make a comment on possible implementation?
    -> In the QG formalism, one can include bottom topography in the RHS of the QG elliptic equation (1b). In terms of implementation, one has to add this term in the 'compute_q_from_psi(self)' method in the 'qgm.py' script.
  
    13) line 114 `smoothness indicator' is too vague. The reader needs a little more explanation
    -> We will provide additional details on this in the revised version.
  
    15) What is the point of the last sentence at the bottom of page 10?
    -> The purpose of this sentence is to inform the reader that this code benefits from automatic differentiation. We do not demonstrate this feature in the present code. We will revise this sentence to make it clearer.
  
  Citation: https://doi.org/10.5194/egusphere-2023-1715-AC1
RC2:
'Comment on egusphere-2023-1715', Anonymous Referee #2, 09 Oct 2023
The authors describe a new numerical solver for the multi-layer quasi-geostrophic equations. The original contributions are of several natures :
numerical : choice of staggering, treatment of immersed boundaries via the capacitance method, WENO advection of potential vorticity without explicit dissipation

computational : concise, GPU-enabled implementation in a high-level language

The paper is very clear and well-written. However it would benefit from a more in-depth discussion of some merits and limitations of the work. Especially :
Regarding the capacitance method, what is its its current and expected numerical cost ? If I am correct one solves a dense linear problem with N unknowns, N being the number of fluid-solid cell edges. Does this cost ∼ N^2 ? For a smooth boundary, one expects N ∼ ( Δ x )^ -1 with Δ x the resolution, so if indeed the cost scales like ∼ N^2 ∼ ( Δ x )^-2 this would not exceed the cost of the rest of the scheme as resolution increases. However realistic coastlines (mentioned line 342) are not smooth, so N may grow faster than N ∼ ( Δ x )^-1 . Could this limit the applicability of the method to more realistic situations ?

Similarly, how far into the hierarchy of models would the presented method still work ? It seems it cannot handle curvilinear coordinates / non-constant metric terms due to the FFT-based solver. Can the authors better delineate the limitations of the method ?

How important is the issue of performance ? Line 235 a single CPU core is about 10x slower than an entire GPU. But usually one would compare full nodes. Is this ratio an indicator of good GPU performance, good CPU performance, neither, both ? Is there room for improvement on the GPU, CPU ?

Vortex shear instabilty : the results shown strongly support the idea that the QG equations are solved correctly and accurately, but this remains qualitative. Due to staggering and despite high-order transport, the whole method is expected to be formally second-order accurate, right ? One could compute the growth rate of the unstable linear mode and discuss its accuracy vs resolution ? Maybe the authors can degrade the order of advection to discuss the added value of high-order advection ?

Double gyre : according to the authors, the challenge is to see an eastward jet emerge in the absence of eddy parameterization. However a single experiment is run. Presumably at low enough resolution the eastward jet disappears. The authors could easily run coarser-resolution experiments and mention the minimal resolution required to have an eastward jet. This would be useful for future comparisons with other methods.

Implementation in a high-level language : can conclusions be drawn for more complex modelling systems ? Independently from the specific numerical choices, is the path followed by the authors doable because the QG system is fairly simple, or is it promising also for, say, Boussinesq ocean models with parameterized vertical mixing ?

I recommend a minor revision that takes into account the above questions and remarks to increase the value of the paper for the readership.

Minor remarks
I could run the provided code but there is no README despite the statement in the text. For the sake of reproducibility it would be useful to provide a conda environment file or similar, with specific versions of python/pytorch/numpy etc.
Citation: https://doi.org/10.5194/egusphere-2023-1715-RC2
- AC2: 'Reply on RC2', Louis Thiry, 09 Nov 2023
  
  We would like to express our gratitude to Referee #2 for their detailed and relevant comments and suggestions.
  
  We address the referee's questions below:
  
  1. Indeed, the capacitance matrix method involves solving a dense linear problem with N unknowns, where N is the number of boundary irregular points (I). This cost scales as N^2, with N^3 pre-computations (capacitance matrix inversion). In practice, the major limitation is in terms of memory. On the laptop utilized, we were able to use up to N=15,000, allowing us to run simulations akin to the North Atlantic at a resolution of 4km. At this resolution, the ageostrophic effects become significant, and the QG hypothesis might no longer be relevant.
  
  2. This is indeed another limitation of the solver: it cannot handle curvilinear coordinates or non-constant metric terms. To address this, we would need to employ a multigrid iterative solver. However, the presented method can still provide an initial estimate for the iterative solver using the mean dx and dy as constant metric terms.
  
  3. The computational efficiency of the QG equations makes them highly practical. We emphasize here the good performance of our solver. To facilitate a fairer comparison, we will include a comparison with full nodes and provide the number of cores for each node, as well as their power consumption in Watts. There is room for improvement for CPU parallelization in our code, as the native parallelism in PyTorch is not as efficient in our case. We will include this in the paper.
  
  4. Vortex shear instability:
  
      a) Correct, thank you for pointing out this aspect. The two-point perpendicular gradient operator is of order two. However, this operator is applied to the stream function, which is the smoothest field, making the benefits of higher-order schemes less evident. We utilize high-order nonlinear reconstruction on the potential vorticity field, which is nonsmooth (boundary currents, eddies, filaments). We will include this discussion in the revision.
  
      b) This is a valuable suggestion, and we will include it in the revised version.
  
      c) This is another insightful suggestion. We will degrade the solution from order 5 to 3 to compare and demonstrate the benefits of high-order advection. Note that for hardware considerations (memory or computation bottlenecks), fifth-order and third-order WENO schemes take nearly the same computational time (see Roullet and Gaillard 2022, JAMES).
  
  5. Double gyre: We will conduct experiments at lower resolutions (27km, 40km, and 53km), including plots to demonstrate how the eastward jet is sensitive to resolution, facilitating comparisons with other methods.
  
  6. We are of the opinion that more complex modeling systems can be implemented in high-level languages like Python. Our plan is to extend the presented approach to shallow-water equations and subsequently to primitive equations. The major advantage is the seamless parallelism offered by GPUs, enabling us to write code that closely aligns with the continuous equations.
  
  Citation: https://doi.org/10.5194/egusphere-2023-1715-AC2

Peer review completion

AR: Author's response | RR: Referee report | ED: Editor decision | EF: Editorial file upload

AR by Louis Thiry on behalf of the Authors (07 Dec 2023) Author's response Author's tracked changes Manuscript

ED: Referee Nomination & Report Request started (18 Dec 2023) by Deepak Subramani

RR by Anonymous Referee #1 (18 Dec 2023)

ED: Publish as is (06 Jan 2024) by Deepak Subramani

AR by Louis Thiry on behalf of the Authors (10 Jan 2024) Manuscript

Short summary

We present a new way of solving the quasi-geostrophic (QG) equations, a simple set of equations describing ocean dynamics. Our method is solely based on the numerical methods used to solve the equations and requires no parameter tuning. Moreover, it can handle non-rectangular geometries, opening the way to study QG equations on realistic domains. We release a PyTorch implementation to ease future machine-learning developments on top of the presented method.