the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 27 Sep 2023
MQGeometry-1.0: a multi-layer quasi-geostrophic solver on non-rectangular geometries
Abstract. This paper presents MQGeometry, a multi-layer quasi-geostrophic (QG) equations solver for non-rectangul 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 stream-function (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 upwind-biased 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 non-rectangular geometries using the capacitance matrix method. We validate our solver on a vortex-shear instability test case in a circular domain, a vortex-wall interaction test-case, and on an idealized wind-driven double-gyre configuration in a octogonal domain at a eddy-permitting resolution. We release a concise, efficient, and auto-differentiable PyTorch implementation of our method to facilitate future developments upon this new discretization, e.g. machine learning parameterization or data-assimilation techniques.
-
Notice on discussion status
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
-
Preprint
(2461 KB)
-
The requested preprint has a corresponding peer-reviewed 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 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
-
AC1: 'Reply on RC1', Louis Thiry, 09 Nov 2023
-
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
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
-
AC1: 'Reply on RC1', Louis Thiry, 09 Nov 2023
-
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
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
Cited
1 citations as recorded by crossref.
Louis Thiry
Guillaume Roullet
Etienne Mémin
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(2461 KB) - Metadata XML