the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
MQGeometry1.0: a multilayer quasigeostrophic solver on nonrectangular geometries
Abstract. This paper presents MQGeometry, a multilayer quasigeostrophic (QG) equations solver for nonrectangul ar geometries. We advect the potential voriticity (PV) with finite volumes to ensure global PV conservation thanks to a staggered discretization of the PV and streamfunction (SF). Thanks to this staggering, the PV is defined inside the domain, removing the need for defining the PV on the domain's boundary. We compute PV fluxes with upwindbiased interpolations whose implicit dissipation replaces the usual explicit (hyper)viscous dissipation. The presented discretization does not require the tuning of any additional parameter, e.g. additional eddy viscosity. We solve the QG elliptic equation with a fast discrete sine transform spectral solver on rectangular geometry. We extend this fast solver to nonrectangular geometries using the capacitance matrix method. We validate our solver on a vortexshear instability test case in a circular domain, a vortexwall interaction testcase, and on an idealized winddriven doublegyre configuration in a octogonal domain at a eddypermitting resolution. We release a concise, efficient, and autodifferentiable PyTorch implementation of our method to facilitate future developments upon this new discretization, e.g. machine learning parameterization or dataassimilation techniques.

Notice on discussion status
The requested preprint has a corresponding peerreviewed final revised paper. You are encouraged to refer to the final revised version.

Preprint
(2461 KB)

The requested preprint has a corresponding peerreviewed final revised paper. You are encouraged to refer to the final revised version.
 Preprint
(2461 KB)  Metadata XML
 BibTeX
 EndNote
 Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed

RC1: 'Comment on egusphere20231715', Anonymous Referee #1, 06 Oct 2023
The paper presents a new implementation of a numerical solver for the multilayer QG equations in nonrectangular 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 wellwritten 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 adhoc...': 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 9395 I cannot understand the sentence `Moreover...'
13) line 114 `smoothness indicator' is too vague. The reader needs a little more explanation
14) l 118119: Please keep the same word order: TVD RungeKutta or RungueKutta 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: fplane > $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/egusphere20231715RC1 
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/egusphere20231715AC1

AC1: 'Reply on RC1', Louis Thiry, 09 Nov 2023

RC2: 'Comment on egusphere20231715', Anonymous Referee #2, 09 Oct 2023
The authors describe a new numerical solver for the multilayer quasigeostrophic 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, GPUenabled implementation in a highlevel language
The paper is very clear and wellwritten. However it would benefit from a more indepth 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 fluidsolid 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 / nonconstant metric terms due to the FFTbased 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 highorder transport, the whole method is expected to be formally secondorder 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 highorder 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 coarserresolution 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 highlevel 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/egusphere20231715RC2 
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 precomputations (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 nonconstant 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 twopoint 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 higherorder schemes less evident. We utilize highorder 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 highorder advection. Note that for hardware considerations (memory or computation bottlenecks), fifthorder and thirdorder 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 highlevel languages like Python. Our plan is to extend the presented approach to shallowwater 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/egusphere20231715AC2
Interactive discussion
Status: closed

RC1: 'Comment on egusphere20231715', Anonymous Referee #1, 06 Oct 2023
The paper presents a new implementation of a numerical solver for the multilayer QG equations in nonrectangular 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 wellwritten 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 adhoc...': 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 9395 I cannot understand the sentence `Moreover...'
13) line 114 `smoothness indicator' is too vague. The reader needs a little more explanation
14) l 118119: Please keep the same word order: TVD RungeKutta or RungueKutta 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: fplane > $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/egusphere20231715RC1 
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/egusphere20231715AC1

AC1: 'Reply on RC1', Louis Thiry, 09 Nov 2023

RC2: 'Comment on egusphere20231715', Anonymous Referee #2, 09 Oct 2023
The authors describe a new numerical solver for the multilayer quasigeostrophic 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, GPUenabled implementation in a highlevel language
The paper is very clear and wellwritten. However it would benefit from a more indepth 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 fluidsolid 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 / nonconstant metric terms due to the FFTbased 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 highorder transport, the whole method is expected to be formally secondorder 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 highorder 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 coarserresolution 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 highlevel 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/egusphere20231715RC2 
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 precomputations (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 nonconstant 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 twopoint 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 higherorder schemes less evident. We utilize highorder 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 highorder advection. Note that for hardware considerations (memory or computation bottlenecks), fifthorder and thirdorder 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 highlevel languages like Python. Our plan is to extend the presented approach to shallowwater 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/egusphere20231715AC2
Peer review completion
Journal article(s) based on this preprint
Viewed
HTML  XML  Total  BibTeX  EndNote  

197  85  17  299  13  7 
 HTML: 197
 PDF: 85
 XML: 17
 Total: 299
 BibTeX: 13
 EndNote: 7
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1
Louis Thiry
Guillaume Roullet
Etienne Mémin
The requested preprint has a corresponding peerreviewed final revised paper. You are encouraged to refer to the final revised version.
 Preprint
(2461 KB)  Metadata XML