the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
The neXtSIM-DG dynamical core: A Framework for Higher-order Finite Element Sea Ice Modeling
Véronique Dansereau
Christian Lessig
Piotr Minakowski
Abstract. The ability of numerical sea ice models to reproduce localized deformation features associated with fracture processes is key for an accurate representation of the ice dynamics and of dynamically coupled physical processes in the Arctic and Antarctic. Equally key is the capacity of these models to minimize the numerical diffusion stemming from the advection of these features, to ensure that the associated strong gradients persist in time, without the need to unphysically re-inject energy for re-localization. To control diffusion and improve the approximation quality, we present a new numerical core for the dynamics of sea ice that is based on higher order finite element discretizations for the momentum equation and higher order discontinuous Galerkin methods for the advection. The mathematical properties of this core are discussed and a detailed description of an efficient shared memory parallel implementation is given. In addition, we present different numerical tests and apply the new framework to a benchmark problem to quantify the advantages of the higher order discretization. These tests are based on Hibler’s viscous-plastic sea ice model, but the implementation of the developed framework in the context of other physical models reproducing a strong localization of the deformation are possible.
-
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
(15385 KB)
-
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(15385 KB) - BibTeX
- EndNote
- Final revised paper
Journal article(s) based on this preprint
Thomas Richter et al.
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2023-391', Sergey Danilov, 28 Mar 2023
It is a well written paper documenting a finite-element based approach to modeling sea ice dynamics in the framework of viscous plastic rheology. The approach is based on quadrilateral elements and explores both bilinear and bi-quadratic representation for sea ice velocity. The tracers are represented using discontinuous elements, but the authors find that the sensitivity to the accuracy of tracer advection is not very high. I recommend publishing this manuscript after minor revision.
I would recommend to briefly discuss the case of nonconforming (CR) element. According to earlier studies, it provides a very high resolution, and will similarly perform very close to the bi-quadratic case. The computational load is in this case seemingly lower. Also, I would suggest to add numbers allowing one to judge about the computational time required by bilinear and bi-quadratic cases. The first one largely corresponds to the approach taken by CICE, so the numbers will be helpful to judge on how the common type discretization is related to the higher-order one.
Minor points:
line 14 'It is, for example, an important part ....' -- It contributes importantly to the global ...15 circulation
62 'to limit diffusion'? Incremental remapping helps to ensure positivity. It relies on limiting, so it is not clear to what an extent it is limiting diffusion.
83 It is the mean height (volume per unit area)
117 remove b
120 It is worthwhile to mention that this is never achieved in practice and mEVP commonly deals with non-converged solutions.
Section 3.1. The discussion of limitations due to advection is not really relevant. First, \Delta t is generally governed by the ocean model, where velocities are larger, so sea ice CFL will not be an issue. Second, there are internal stresses that depend on the mean thickness and concentration, and this may lead to a limitation
(for coupled equations) that is more demanding.130 Sea ice was sometimes run with larger time step than ocean, but I am not aware about the situation described here.
133 Why this range? It will obviously depend on mesh resolution.
135 This is the limitations due to advection, it may matter for high-order in principle, but as I wrote, it is commonly ocean that matters most, and there is a wave-type limitation due to plastic response.
140 To be consistent ... -- MPAS, ICON and FESOM communities occupy a substantial part of climate modeling, and they rely on different meshes
141 unstructured --> distorted
144 for a depiction --> for an illustration
Expressions (12) and (13): the lower indices start from 0 in (13), and 1 in (12)
186 Which order of RK is used?
212 The Babushka-Brezzi condition is a subject of many publications and is well know, so at least include (see, e.g. Ern ...)
Remark 1 -- Why 'Optimality...'? Say before that the discretization on quads involves a lot of DoFs and computations, and a question arises on its optimality with respect to triangular meshes where functional spaces seem smaller.
256 by Gaussian325 a triangle -- an element
Table 1 Ice concentration?
450 is strongly affected ? I do not think it is the case. As soon as stability is ensured, sensitivity is generally very moderate.
466 THis
5.2.3. I do not think this should go in the main text, move it to an appendix
Figure 14, caption '10 momentum steps' --- 10 or 100?
5.3.1. The efficiency on shared memory level is perhaps not surprising given the measures described. The real challenge, however, is the MPI parallelism, because of the large number of substeps in mEVP.
Sergey Danilov
Citation: https://doi.org/10.5194/egusphere-2023-391-RC1 -
AC1: 'Reply on RC1', Thomas Richter, 05 Apr 2023
Dear Sergey,
many thanks for your comments that we will all gladly take into account. But for the further discussion I would already like to post some results regarding the computational times:
dG(0) dG(1) dG(2) cG(1) 47 s 51 s 60 s cG(2) 89 s 91 s 102 s This is the benchmark setup on [0, 2 days], a 4km mesh, 120s time-step size, 100 mEVP iterations using 32 Cores on an AMD Epyc CPU. Times include everything apart from i/o. About 90% (more for dG(0)) is the attributed to the momentum equation, the remainder to the advection.
Due to increased local work, cG(2) is more efficient in terms of vectorization and parallelization. This explains the rather good scaling, less than twice the comp. times compared to cG(1).
Best,
Thomas
Citation: https://doi.org/10.5194/egusphere-2023-391-AC1
-
AC1: 'Reply on RC1', Thomas Richter, 05 Apr 2023
-
RC2: 'Comment on egusphere-2023-391', Anonymous Referee #2, 10 Apr 2023
The authors present a finite-element model of sea-ice under elastic-visco-plastic approximation. Specifically, they chose a dG advection scheme for sea-ice concentration, and a mixed approach for the momentum equation, with continuous Galerkin method for all terms except the stress tensor, which was discretized with dG method. In the numerical tests, the authors focus on how different order cG and dG discretizations affect the sea-ice shear deformation, and study computational performance of their code parallelized using shared memory paradigm. Despite the increased computational cost, higher order discretization seem to provide superior results and much smaller diffusion. I am glad that authors consider moving to even higher basis polynomial orders, as the increased computational per-element intensity can lead to better parallel scalability results, offsetting the increased cost.
The paper is well written and contains sufficient math background, numerical tests and physics validation to merit publication. I do not have recommendation for revisions.
Citation: https://doi.org/10.5194/egusphere-2023-391-RC2
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2023-391', Sergey Danilov, 28 Mar 2023
It is a well written paper documenting a finite-element based approach to modeling sea ice dynamics in the framework of viscous plastic rheology. The approach is based on quadrilateral elements and explores both bilinear and bi-quadratic representation for sea ice velocity. The tracers are represented using discontinuous elements, but the authors find that the sensitivity to the accuracy of tracer advection is not very high. I recommend publishing this manuscript after minor revision.
I would recommend to briefly discuss the case of nonconforming (CR) element. According to earlier studies, it provides a very high resolution, and will similarly perform very close to the bi-quadratic case. The computational load is in this case seemingly lower. Also, I would suggest to add numbers allowing one to judge about the computational time required by bilinear and bi-quadratic cases. The first one largely corresponds to the approach taken by CICE, so the numbers will be helpful to judge on how the common type discretization is related to the higher-order one.
Minor points:
line 14 'It is, for example, an important part ....' -- It contributes importantly to the global ...15 circulation
62 'to limit diffusion'? Incremental remapping helps to ensure positivity. It relies on limiting, so it is not clear to what an extent it is limiting diffusion.
83 It is the mean height (volume per unit area)
117 remove b
120 It is worthwhile to mention that this is never achieved in practice and mEVP commonly deals with non-converged solutions.
Section 3.1. The discussion of limitations due to advection is not really relevant. First, \Delta t is generally governed by the ocean model, where velocities are larger, so sea ice CFL will not be an issue. Second, there are internal stresses that depend on the mean thickness and concentration, and this may lead to a limitation
(for coupled equations) that is more demanding.130 Sea ice was sometimes run with larger time step than ocean, but I am not aware about the situation described here.
133 Why this range? It will obviously depend on mesh resolution.
135 This is the limitations due to advection, it may matter for high-order in principle, but as I wrote, it is commonly ocean that matters most, and there is a wave-type limitation due to plastic response.
140 To be consistent ... -- MPAS, ICON and FESOM communities occupy a substantial part of climate modeling, and they rely on different meshes
141 unstructured --> distorted
144 for a depiction --> for an illustration
Expressions (12) and (13): the lower indices start from 0 in (13), and 1 in (12)
186 Which order of RK is used?
212 The Babushka-Brezzi condition is a subject of many publications and is well know, so at least include (see, e.g. Ern ...)
Remark 1 -- Why 'Optimality...'? Say before that the discretization on quads involves a lot of DoFs and computations, and a question arises on its optimality with respect to triangular meshes where functional spaces seem smaller.
256 by Gaussian325 a triangle -- an element
Table 1 Ice concentration?
450 is strongly affected ? I do not think it is the case. As soon as stability is ensured, sensitivity is generally very moderate.
466 THis
5.2.3. I do not think this should go in the main text, move it to an appendix
Figure 14, caption '10 momentum steps' --- 10 or 100?
5.3.1. The efficiency on shared memory level is perhaps not surprising given the measures described. The real challenge, however, is the MPI parallelism, because of the large number of substeps in mEVP.
Sergey Danilov
Citation: https://doi.org/10.5194/egusphere-2023-391-RC1 -
AC1: 'Reply on RC1', Thomas Richter, 05 Apr 2023
Dear Sergey,
many thanks for your comments that we will all gladly take into account. But for the further discussion I would already like to post some results regarding the computational times:
dG(0) dG(1) dG(2) cG(1) 47 s 51 s 60 s cG(2) 89 s 91 s 102 s This is the benchmark setup on [0, 2 days], a 4km mesh, 120s time-step size, 100 mEVP iterations using 32 Cores on an AMD Epyc CPU. Times include everything apart from i/o. About 90% (more for dG(0)) is the attributed to the momentum equation, the remainder to the advection.
Due to increased local work, cG(2) is more efficient in terms of vectorization and parallelization. This explains the rather good scaling, less than twice the comp. times compared to cG(1).
Best,
Thomas
Citation: https://doi.org/10.5194/egusphere-2023-391-AC1
-
AC1: 'Reply on RC1', Thomas Richter, 05 Apr 2023
-
RC2: 'Comment on egusphere-2023-391', Anonymous Referee #2, 10 Apr 2023
The authors present a finite-element model of sea-ice under elastic-visco-plastic approximation. Specifically, they chose a dG advection scheme for sea-ice concentration, and a mixed approach for the momentum equation, with continuous Galerkin method for all terms except the stress tensor, which was discretized with dG method. In the numerical tests, the authors focus on how different order cG and dG discretizations affect the sea-ice shear deformation, and study computational performance of their code parallelized using shared memory paradigm. Despite the increased computational cost, higher order discretization seem to provide superior results and much smaller diffusion. I am glad that authors consider moving to even higher basis polynomial orders, as the increased computational per-element intensity can lead to better parallel scalability results, offsetting the increased cost.
The paper is well written and contains sufficient math background, numerical tests and physics validation to merit publication. I do not have recommendation for revisions.
Citation: https://doi.org/10.5194/egusphere-2023-391-RC2
Peer review completion
Journal article(s) based on this preprint
Thomas Richter et al.
Thomas Richter et al.
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 330 | 110 | 15 | 455 | 7 | 3 |
- HTML: 330
- PDF: 110
- XML: 15
- Total: 455
- BibTeX: 7
- EndNote: 3
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(15385 KB) - Metadata XML