-
Local Time-Stepping for the Shallow Water Equations using CFL Optimized Forward-Backward Runge-Kutta Schemes
Authors:
Jeremy R. Lilly,
Giacomo Capodaglio,
Darren Engwirda,
Robert L. Higdon,
Mark R. Petersen
Abstract:
The Courant-Friedrichs-Lewy (CFL) condition is a well known, necessary condition for the stability of explicit time-stepping schemes that effectively places a limit on the size of the largest admittable time-step for a given problem. We formulate and present a new local time-stepping (LTS) scheme optimized, in the CFL sense, for the shallow water equations (SWEs). This new scheme, called FB-LTS, i…
▽ More
The Courant-Friedrichs-Lewy (CFL) condition is a well known, necessary condition for the stability of explicit time-stepping schemes that effectively places a limit on the size of the largest admittable time-step for a given problem. We formulate and present a new local time-stepping (LTS) scheme optimized, in the CFL sense, for the shallow water equations (SWEs). This new scheme, called FB-LTS, is based on the CFL optimized forward-backward Runge-Kutta schemes from Lilly et al. (2023). We show that FB-LTS maintains exact conservation of mass and absolute vorticity when applied to the TRiSK spatial discretization (Ringler et al., 2010), and provide numerical experiments showing that it retains the temporal order of the scheme on which it is based (second order). In terms of computational performance, we show that when applied to a real-world test case on a highly-variable resolution mesh, the MPAS-Ocean implementation of FB-LTS is up to 10 times faster than the classical four-stage, fourth-order Runge-Kutta method (RK4), and 2.3 times faster than an existing strong stability preserving Runge-Kutta based LTS scheme (LTS3). Despite this significant increase in efficiency, the solutions produced by FB-LTS are qualitatively equivalent to those produced by both RK4 and LTS3.
△ Less
Submitted 16 May, 2024;
originally announced May 2024.
-
CFL Optimized Forward-Backward Runge-Kutta Schemes for the Shallow Water Equations
Authors:
Jeremy R. Lilly,
Darren Engwirda,
Giacomo Capodaglio,
Robert L. Higdon,
Mark R. Petersen
Abstract:
We present the formulation and optimization of a Runge-Kutta-type time-stepping scheme for solving the shallow water equations, aimed at substantially increasing the effective allowable time-step over that of comparable methods. This scheme, called FB-RK(3,2), uses weighted forward-backward averaging of thickness data to advance the momentum equation. The weights for this averaging are chosen with…
▽ More
We present the formulation and optimization of a Runge-Kutta-type time-stepping scheme for solving the shallow water equations, aimed at substantially increasing the effective allowable time-step over that of comparable methods. This scheme, called FB-RK(3,2), uses weighted forward-backward averaging of thickness data to advance the momentum equation. The weights for this averaging are chosen with an optimization process that employs a von Neumann-type analysis, ensuring that the weights maximize the admittable Courant number. Through a simplified local truncation error analysis and numerical experiments, we show that the method is at least second order in time for any choice of weights and exhibits low dispersion and dissipation errors for well-resolved waves. Further, we show that an optimized FB-RK(3,2) can take time-steps up to 2.8 times as large as a popular three-stage, third-order strong stability preserving Runge-Kutta method in a quasi-linear test case. In fully nonlinear shallow water test cases relevant to oceanic and atmospheric flows, FB-RK(3,2) outperforms SSPRK3 in admittable time-step by factors between roughly between 1.6 and 2.2, making the scheme approximately twice as computationally efficient with little to no effect on solution quality.
△ Less
Submitted 20 June, 2023;
originally announced June 2023.
-
A scalable domain decomposition method for FEM discretizations of nonlocal equations of integrable and fractional type
Authors:
Manuel Klar,
Giacomo Capodaglio,
Marta D'Elia,
Christian Glusa,
Max Gunzburger,
Christian Vollmann
Abstract:
Nonlocal models allow for the description of phenomena which cannot be captured by classical partial differential equations. The availability of efficient solvers is one of the main concerns for the use of nonlocal models in real world engineering applications. We present a domain decomposition solver that is inspired by substructuring methods for classical local equations. In numerical experiment…
▽ More
Nonlocal models allow for the description of phenomena which cannot be captured by classical partial differential equations. The availability of efficient solvers is one of the main concerns for the use of nonlocal models in real world engineering applications. We present a domain decomposition solver that is inspired by substructuring methods for classical local equations. In numerical experiments involving finite element discretizations of scalar and vectorial nonlocal equations of integrable and fractional type, we observe improvements in solution time of up to 14.6x compared to commonly used solver strategies.
△ Less
Submitted 31 May, 2023;
originally announced June 2023.
-
Simulating sea-ice deformation in viscous-plastic sea-ice models with CD-grids
Authors:
C. Mehlmann,
G. Capodaglio,
S. Danilov
Abstract:
Linear Kinematic Features (LKFs) are found everywhere in the Arctic sea-ice cover. They are strongly localized deformations often associated with the formation of leads and pressure ridges. Viscous-plastic sea-ice models start to produce LKFs at high spatial grid resolution, typically with a grid spacing below 5 km. A recent study showed that the placement of the variables on the grid plays an imp…
▽ More
Linear Kinematic Features (LKFs) are found everywhere in the Arctic sea-ice cover. They are strongly localized deformations often associated with the formation of leads and pressure ridges. Viscous-plastic sea-ice models start to produce LKFs at high spatial grid resolution, typically with a grid spacing below 5 km. A recent study showed that the placement of the variables on the grid plays an important role for the number of simulated LKFs. The study found that a nonconforming finite element discretization with a CD-grid placement (CD1) resolves more LKFs per degree of freedom compared to more common A,B and C-grids. A new CD-grid formulation (CD2) has just been proposed based on a conforming subgrid discretization. To analyze the resolution properties of the new CD2 approach we evaluate runs from different models (e.g FESOM, MPAS) on a benchmark problem using quadrilateral, hexagonal and triangular meshes. We found that the CD1 setup simulates more deformation structure than the CD2 approximation. This highlights the importance of the type of spatial discretization for the simulation of LKFs. Due to the higher number of degrees of freedom both CD-grids resolve more LKFs than traditional A,B and C-grids at fixed mesh level. This is an appealing feature as high spatial mesh resolution is needed in viscous-plastic sea-ice models to simulate LKFs.
△ Less
Submitted 17 April, 2023;
originally announced April 2023.
-
Storm Surge Modeling as an Application of Local Time-stepping in MPAS-Ocean
Authors:
Jeremy R. Lilly,
Giacomo Capodaglio,
Mark R. Petersen,
Steven R. Brus,
Darren Engwirda,
Robert L. Higdon
Abstract:
This paper presents the first scientific application of local time-stepping (LTS) schemes in the Model for Prediction Across Scales-Ocean (MPAS-O). We use LTS schemes in a single-layer, global ocean model that predicts the storm surge around the eastern coast of the United States during Hurricane Sandy. The variable-resolution meshes used are of unprecedentedly high resolution in MPAS-O, containin…
▽ More
This paper presents the first scientific application of local time-stepping (LTS) schemes in the Model for Prediction Across Scales-Ocean (MPAS-O). We use LTS schemes in a single-layer, global ocean model that predicts the storm surge around the eastern coast of the United States during Hurricane Sandy. The variable-resolution meshes used are of unprecedentedly high resolution in MPAS-O, containing cells as small as 125 meters wide in Delaware Bay. It is shown that a particular, third-order LTS scheme (LTS3) produces sea-surface height (SSH) solutions that are of comparable quality to solutions produced by the classical four-stage, fourth-order Runge-Kutta method (RK4) with a uniform time step on the same meshes. Furthermore, LTS3 is up to 35% faster in the best cases, showing that LTS schemes are viable for use in MPAS-O with the added benefit of substantially less computational cost. The results of these performance experiments inform us of the requirements for efficient mesh design for LTS schemes. In particular, we see that for LTS to be efficient on a given mesh, it is important to have enough cells using the coarse time-step relative to those using the fine time-step, typically at least 1:5.
△ Less
Submitted 27 July, 2022;
originally announced July 2022.
-
An asymptotically compatible coupling formulation for nonlocal interface problems with jumps
Authors:
Christian Glusa,
Marta D'Elia,
Giacomo Capodaglio,
Max Gunzburger,
Pavel B. Bochev
Abstract:
We introduce a mathematically rigorous formulation for a nonlocal interface problem with jumps and propose an asymptotically compatible finite element discretization for the weak form of the interface problem. After proving the well-posedness of the weak form, we demonstrate that solutions to the nonlocal interface problem converge to the corresponding local counterpart when the nonlocal data are…
▽ More
We introduce a mathematically rigorous formulation for a nonlocal interface problem with jumps and propose an asymptotically compatible finite element discretization for the weak form of the interface problem. After proving the well-posedness of the weak form, we demonstrate that solutions to the nonlocal interface problem converge to the corresponding local counterpart when the nonlocal data are appropriately prescribed. Several numerical tests in one and two dimensions show the applicability of our technique, its numerical convergence to exact nonlocal solutions, its convergence to the local limit when the horizons vanish, and its robustness with respect to the patch test.
△ Less
Submitted 14 March, 2022;
originally announced March 2022.
-
An unstructured CD-grid variational formulation for sea ice dynamics
Authors:
Giacomo Capodaglio,
Mark R. Petersen,
Adrian K. Turner,
Andrew F. Roberts
Abstract:
For the numerical simulation of earth system models, Arakawa grids are largely employed. A quadrilateral mesh is assumed for their classical definition, and different types of grids are identified depending on the location of the discretized quantities. The B-grid has both velocity components at the center of a cell, the C-grid places the velocity components on the edges in a staggered fashion, an…
▽ More
For the numerical simulation of earth system models, Arakawa grids are largely employed. A quadrilateral mesh is assumed for their classical definition, and different types of grids are identified depending on the location of the discretized quantities. The B-grid has both velocity components at the center of a cell, the C-grid places the velocity components on the edges in a staggered fashion, and the D-grid is a ninety-degree rotation of a C-grid. Historically, B-grid formulations of sea ice dynamics have been dominant because they have matched the grid type used by ocean models. In recent years, as ocean models have increasingly progressed to C-grids, sea ice models have followed suit on quadrilateral meshes, but few if any implementations of unstructured C-grid sea ice models have been developed. In this work, we present an unstructured CD-grid type formulation of the elastic-viscous-plastic rheology, where the velocity unknowns are located at the edges, rather than at the vertices, as in the B-grid. The notion of a CD-grid has been recently introduced and assumes that the velocity components are co-located at the edges. The mesh cells in our analysis have $n$ sides, with $n$ greater than or equal to four. Numerical results are included to investigate the features of the proposed method. Our framework of choice is the Model for Prediction Across Scales (MPAS) within E3SM, the climate model of the U.S. Department of Energy, although our approach is general and could be applied to other models as well. While MPAS-Seaice is currently defined on a B-grid, MPAS-Ocean runs on a C-grid, hence interpolation operators are heavily used when coupled simulations are performed. The discretization introduced here aims at transitioning the dynamics of MPAS-Seaice to a CD-grid, to ultimately facilitate improved coupling with MPAS-Ocean and reduce numerical errors associated with this communication.
△ Less
Submitted 6 July, 2022; v1 submitted 30 December, 2021;
originally announced December 2021.
-
Local time stepping for the shallow water equations in MPAS
Authors:
Giacomo Capodaglio,
Mark Petersen
Abstract:
We assess the performance of a set of local time-stepping (LTS) schemes for the shallow water equations implemented in the Model for Prediction Across Scales (MPAS). The goal of LTS is to speed up the simulation by allowing different time-steps on different regions of the computational grid. The LTS schemes considered here were originally introduced by Hoang et al. (J. Comput. Phys., Vol. 382, p.1…
▽ More
We assess the performance of a set of local time-stepping (LTS) schemes for the shallow water equations implemented in the Model for Prediction Across Scales (MPAS). The goal of LTS is to speed up the simulation by allowing different time-steps on different regions of the computational grid. The LTS schemes considered here were originally introduced by Hoang et al. (J. Comput. Phys., Vol. 382, p.152-176, 2019), who laid out the mathematical foundation of the methods. Here, the authors take on the task of presenting a fast, efficient and scalable parallel implementation of these LTS methods on high performance computing machines, with the aim to provide a recipe for other climate modeling groups that may be interested in employing LTS algorithms in their codes. As a matter of fact, even if MPAS is our framework of choice, our approach is general enough and could be of interest to other groups beyond the MPAS community. Due to their nature, LTS methods possess an inherent load imbalance that needs to be carefully addressed in order to obtain efficient scalability. Even more important is the far from trivial task of computing the right-hand side terms only on specific LTS regions during the time-stepping procedure. An inefficient handling of this task causes a drastic decay of the CPU time performance, making the LTS algorithms practically of no use. The emphasis of the present work is therefore on the computational and parallel aspects of the LTS methods, whose proper handling is crucial to make the methods run faster against existing strategies, such as for instance high-order explicit global time-stepping schemes. This is in fact the ultimate goal of using an LTS procedure and it is the one to which we direct all our optimization efforts.
△ Less
Submitted 22 September, 2021; v1 submitted 14 June, 2021;
originally announced June 2021.
-
Efficient quadrature rules for finite element discretizations of nonlocal equations
Authors:
Eugenio Aulisa,
Giacomo Capodaglio,
Andrea Chierici,
Marta D'Elia
Abstract:
In this paper we design efficient quadrature rules for finite element discretizations of nonlocal diffusion problems with compactly supported kernel functions. Two of the main challenges in nonlocal modeling and simulations are the prohibitive computational cost and the nontrivial implementation of discretization schemes, especially in three-dimensional settings. In this work we circumvent both ch…
▽ More
In this paper we design efficient quadrature rules for finite element discretizations of nonlocal diffusion problems with compactly supported kernel functions. Two of the main challenges in nonlocal modeling and simulations are the prohibitive computational cost and the nontrivial implementation of discretization schemes, especially in three-dimensional settings. In this work we circumvent both challenges by introducing a parametrized mollifying function that improves the regularity of the integrand, utilizing an adaptive integration technique, and exploiting parallelization. We first show that the "mollified" solution converges to the exact one as the mollifying parameter vanishes, then we illustrate the consistency and accuracy of the proposed method on several two- and three-dimensional test cases. Furthermore, we demonstrate the good scaling properties of the parallel implementation of the adaptive algorithm and we compare the proposed method with recently developed techniques for efficient finite element assembly.
△ Less
Submitted 21 January, 2021;
originally announced January 2021.
-
A general framework for substructuring-based domain decomposition methods for models having nonlocal interactions
Authors:
Giacomo Capodaglio,
Marta D'Elia,
Max Gunzburger,
Pavel Bochev,
Manuel Klar,
Christian Vollmann
Abstract:
A rigorous mathematical framework is provided for a substructuring-based domain-decomposition approach for nonlocal problems that feature interactions between points separated by a finite distance. Here, by substructuring it is meant that a traditional geometric configuration for local partial differential equation problems is used in which a computational domain is subdivided into non-overlapping…
▽ More
A rigorous mathematical framework is provided for a substructuring-based domain-decomposition approach for nonlocal problems that feature interactions between points separated by a finite distance. Here, by substructuring it is meant that a traditional geometric configuration for local partial differential equation problems is used in which a computational domain is subdivided into non-overlapping subdomains. In the nonlocal setting, this approach is substructuring-based in the sense that those subdomains interact with neighboring domains over interface regions having finite volume, in contrast to the local PDE setting in which interfaces are lower dimensional manifolds separating abutting subdomains. Key results include the equivalence between the global, single-domain nonlocal problem and its multi-domain reformulation, both at the continuous and discrete levels. These results provide the rigorous foundation necessary for the development of efficient solution strategies for nonlocal domain-decomposition methods.
△ Less
Submitted 26 August, 2020;
originally announced August 2020.
-
An energy-based coupling approach to nonlocal interface problems
Authors:
Giacomo Capodaglio,
Marta D'Elia,
Pavel Bochev,
Max Gunzburger
Abstract:
Nonlocal models provide accurate representations of physical phenomena ranging from fracture mechanics to complex subsurface flows, where traditional partial differential equations fail to capture effects caused by long-range forces at the microscale and mesoscale. However, the application of nonlocal models to problems involving interfaces such as multimaterial simulations and fluid-structure int…
▽ More
Nonlocal models provide accurate representations of physical phenomena ranging from fracture mechanics to complex subsurface flows, where traditional partial differential equations fail to capture effects caused by long-range forces at the microscale and mesoscale. However, the application of nonlocal models to problems involving interfaces such as multimaterial simulations and fluid-structure interaction, is hampered by the lack of a rigorous nonlocal interface theory needed to support numerical developments. In this paper, we use an energy-based approach to develop a mathematically rigorous nonlocal interface theory which provides a physically consistent extension of the classical perfect interface PDE formulation. Numerical examples validate the proposed framework and demonstrate the scope of our theory.
△ Less
Submitted 7 May, 2020; v1 submitted 10 January, 2020;
originally announced January 2020.
-
A computational study of preconditioning techniques for the stochastic diffusion equation with lognormal coefficient
Authors:
Eugenio Aulisa,
Giacomo Capodaglio,
Guoyi Ke
Abstract:
We present a computational study of several preconditioning techniques for the GMRES algorithm applied to the stochastic diffusion equation with a lognormal coefficient discretized with the stochastic Galerkin method. The clear block structure of the system matrix arising from this type of discretization motivates the analysis of preconditioners designed according to a field-splitting strategy of…
▽ More
We present a computational study of several preconditioning techniques for the GMRES algorithm applied to the stochastic diffusion equation with a lognormal coefficient discretized with the stochastic Galerkin method. The clear block structure of the system matrix arising from this type of discretization motivates the analysis of preconditioners designed according to a field-splitting strategy of the stochastic variables. This approach is inspired by a similar procedure used within the framework of physics based preconditioners for deterministic problems, and its application to stochastic PDEs represents the main novelty of this work. Our numerical investigation highlights the superior properties of the field-split type preconditioners over other existing strategies in terms of computational time and stochastic parameter dependence.
△ Less
Submitted 24 October, 2019;
originally announced October 2019.
-
Piecewise polynomial approximation of probability density functions with application to uncertainty quantification for stochastic PDEs
Authors:
Giacomo Capodaglio,
Max Gunzburger
Abstract:
The probability density function (PDF) associated with a given set of samples is approximated by a piecewise-linear polynomial constructed with respect to a binning of the sample space. The kernel functions are a compactly supported basis for the space of such polynomials, i.e. finite element hat functions, that are centered at the bin nodes rather than at the samples, as is the case for the stand…
▽ More
The probability density function (PDF) associated with a given set of samples is approximated by a piecewise-linear polynomial constructed with respect to a binning of the sample space. The kernel functions are a compactly supported basis for the space of such polynomials, i.e. finite element hat functions, that are centered at the bin nodes rather than at the samples, as is the case for the standard kernel density estimation approach. This feature naturally provides an approximation that is scalable with respect to the sample size. On the other hand, unlike other strategies that use a finite element approach, the proposed approximation does not require the solution of a linear system. In addition, a simple rule that relates the bin size to the sample size eliminates the need for bandwidth selection procedures. The proposed density estimator has unitary integral, does not require a constraint to enforce positivity, and is consistent. The proposed approach is validated through numerical examples in which samples are drawn from known PDFs. The approach is also used to determine approximations of (unknown) PDFs associated with outputs of interest that depend on the solution of a stochastic partial differential equation.
△ Less
Submitted 5 December, 2019; v1 submitted 26 June, 2019;
originally announced June 2019.
-
Monolithic coupling of implicit material point method with finite element method
Authors:
Eugenio Aulisa,
Giacomo Capodaglio
Abstract:
A monolithic coupling between the material point method (MPM) and the finite element method (FEM) is presented. The MPM formulation described is implicit, and the exchange of information between particles and background grid is minimized. The reduced information transfer from the particles to the grid improves the stability of the method. Once the residual is assembled, the system matrix is obtain…
▽ More
A monolithic coupling between the material point method (MPM) and the finite element method (FEM) is presented. The MPM formulation described is implicit, and the exchange of information between particles and background grid is minimized. The reduced information transfer from the particles to the grid improves the stability of the method. Once the residual is assembled, the system matrix is obtained by means of automatic differentiation. In such a way, no explicit computation is required and the implementation is considerably simplified. When MPM is coupled with FEM, the MPM background grid is attached to the FEM body and the coupling is monolithic. With this strategy, no MPM particle can penetrate a FEM element, and the need for computationally expensive contact search algorithms used by existing coupling procedures is eliminated. The coupled system can be assembled with a single assembly procedure carried out element by element in a FEM fashion. Numerical results are reported to display the performances and advantages of the methods here discussed.
△ Less
Submitted 27 November, 2018;
originally announced November 2018.
-
Approximation of probability density functions for PDEs with random parameters using truncated series expansions
Authors:
Giacomo Capodaglio,
Max Gunzburger,
Henry P. Wynn
Abstract:
The probability density function (PDF) of a random variable associated with the solution of a partial differential equation (PDE) with random parameters is approximated using a truncated series expansion. The random PDE is solved using two stochastic finite element methods, Monte Carlo sampling and the stochastic Galerkin method with global polynomials. The random variable is a functional of the s…
▽ More
The probability density function (PDF) of a random variable associated with the solution of a partial differential equation (PDE) with random parameters is approximated using a truncated series expansion. The random PDE is solved using two stochastic finite element methods, Monte Carlo sampling and the stochastic Galerkin method with global polynomials. The random variable is a functional of the solution of the random PDE, such as the average over the physical domain. The truncated series are obtained considering a finite number of terms in the Gram-Charlier or Edgeworth series expansions. These expansions approximate the PDF of a random variable in terms of another PDF, and involve coefficients that are functions of the known cumulants of the random variable. To the best of our knowledge, their use in the framework of PDEs with random parameters has not yet been explored.
△ Less
Submitted 23 September, 2020; v1 submitted 1 October, 2018;
originally announced October 2018.
-
Construction of h-refined continuous finite element spaces with arbitrary hanging node configurations and applications to multigrid algorithms
Authors:
Eugenio Aulisa,
Giacomo Capodaglio,
Guoyi Ke
Abstract:
We present a novel approach for the construction of basis functions to be employed in selective or adaptive h-refined finite element applications with arbitrary-level hanging node configurations. Our analysis is not restricted to $1$-irregular meshes, as it is usually done in the literature, allowing our results to be applicable to a broader class of local refinement strategies. The proposed metho…
▽ More
We present a novel approach for the construction of basis functions to be employed in selective or adaptive h-refined finite element applications with arbitrary-level hanging node configurations. Our analysis is not restricted to $1$-irregular meshes, as it is usually done in the literature, allowing our results to be applicable to a broader class of local refinement strategies. The proposed method does not require the solution of any linear system to obtain the constraints necessary to enforce continuity of the basis functions and it can be easily implemented. A mathematical analysis is carried out to prove that the proposed basis functions are continuous and linearly independent. Finite element spaces are then defined as the spanning sets of such functions, and the implementation of a multigrid algorithm built on these spaces is discussed. A spectral analysis of the multigrid algorithm highlights superior convergence properties of the proposed method over existing strategies based on a local smoothing procedure. Finally, linear and nonlinear numerical examples are tested to show the robustness and versatility of the multigrid algorithm.
△ Less
Submitted 1 May, 2018; v1 submitted 27 April, 2018;
originally announced April 2018.
-
Convergence estimates for multigrid algorithms with SSC smoothers and applications to overlapping domain decomposition
Authors:
Eugenio Aulisa,
Giorgio Bornia,
Sara Calandrini,
Giacomo Capodaglio
Abstract:
In this paper we study convergence estimates for a multigrid algorithm with smoothers of successive subspace correction (SSC) type, applied to symmetric elliptic PDEs. First, we revisit a general convergence analysis on a class of multigrid algorithms in a fairly general setting, where no regularity assumptions are made on the solution. In this framework, we are able to explicitly highlight the de…
▽ More
In this paper we study convergence estimates for a multigrid algorithm with smoothers of successive subspace correction (SSC) type, applied to symmetric elliptic PDEs. First, we revisit a general convergence analysis on a class of multigrid algorithms in a fairly general setting, where no regularity assumptions are made on the solution. In this framework, we are able to explicitly highlight the dependence of the multigrid error bound on the number of smoothing steps. For the case of no regularity assumptions, this represents a new addition to the existing theory. Then, we analyze successive subspace correction smoothing schemes for a set of uniform and local refinement applications with either nested or non-nested overlapping subdomains. For these applications, we explicitly derive bounds for the multigrid error, and identify sufficient conditions for these bounds to be independent of the number of multigrid levels. For the local refinement applications, finite element grids with arbitrary hanging nodes configurations are considered. The analysis of these smoothing schemes is cast within the far-reaching multiplicative Schwarz framework.
△ Less
Submitted 27 September, 2017;
originally announced September 2017.