Search

This paper complements the empirical justification of the revised scheme in Part I of this work with a mathematical justification leveraging a semi-discrete analysis framework for assessing the splitting error of process coupling methods. The novelty of the framework is that splitting error is distinguished from the process time integration errors, i.e., the errors caused by discrete time integration of individual processes, leading to expressions that are more easily interpreted utilizing existing physical understanding of the processes that the terms represent. This application of this framework to dust life cycle in EAMv1 showcases such an interpretation, using the leading-order splitting error that results from the framework to confirm (i) that the original EAMv1 scheme artificially strengthens the effect of dry removal processes, and (ii) that the revised splitting reduces that artificial strengthening. While the error analysis framework is presented in the context of the dust life cycle in EAMv1, the framework can be broadly leveraged to evaluate process coupling schemes, both in other physical problems and for any number of processes. This framework will be particularly powerful when the various process implementations support a variety of time integration approaches. Whereas traditional local truncation error approaches require separate consideration of each combination of time integration methods, this framework enables evaluation of coupling schemes independent of particular time integration approaches for each process while still allowing for the incorporation of these specific time integration errors if so desired. The framework also explains how the splitting error terms result from (i) the integration of individual processes in isolation from other processes, and (ii) the choices of input state and timestep size for the isolated integration of processes.

An earlier study evaluating the dust life cycle in the Energy Exascale Earth System Model (E3SM) Atmosphere Model version 1 (EAMv1) has revealed that the simulated global mean dust lifetime is substantially shorter when higher vertical resolution is used, primarily due to significant strengthening of dust dry removal in source regions. This paper demonstrates that the sequential splitting of aerosol emissions, dry removal, and turbulent mixing in the model's time integration loop, especially the calculation of dry removal after surface emissions and before turbulent mixing, is the primary reason for the vertical resolution sensitivity reported in that earlier study. Based on this reasoning, we propose a simple revision to the numerical process coupling scheme, which moves the application of the surface emissions to after dry removal and before turbulent mixing. The revised scheme allows newly emitted particles to be transported aloft by turbulence before being removed from the atmosphere, and hence better resembles the dust life cycle in the real world. Sensitivity experiments are conducted and analyzed to evaluate the impact of the revised coupling on the simulated aerosol climatology in EAMv1.

Highly accurate simulation of plasma transport is needed for the successful design and operation of magnetically confined fusion reactors. Unfortunately, the extreme anisotropy present in magnetized plasmas results in thin boundary layers that are expensive to resolve. This work investigates how mesh refinement strategies might reduce that expense to allow for more efficient simulation. It is first verified that higher order discretization only realizes the proper rate of convergence once the mesh resolves the thin boundary layer, motivating the focusing of refinement on the boundary layer. Three mesh refinement strategies are investigated: one that focuses the refinement across the layer by using rectangular elements with a ratio equal to the boundary layer width, one that allows for exponential growth in mesh spacing away from the layer, and one adaptive strategy utilizing the established Zienkiewicz and Zhu error estimator. Across 4 two-dimensional test cases with high anisotropy, the adaptive mesh refinement strategy consistently achieves the same accuracy as uniform refinement using orders of magnitude less degrees of freedom. In the test case where the magnetic field is aligned with the mesh, the other refinement strategies also show substantial improvement in efficiency. This work also includes a discussion generalizing the results to larger magnetic anisotropy ratios and to three-dimensional problems. It is shown that isotropic mesh refinement requires degrees of freedom on the order of either the layer width (2D) or the square of the layer width (3D), whereas anisotropic refinement requires a number on the order of the log of layer width for all dimensions. It is also shown that the number of conjugate gradient iterations scales as a power of layer width when preconditioned with algebraic multigrid, whereas the number is independent of layer width when preconditioned with ILU.

Moving scientific computation from high-performance computing (HPC) and cloud computing (CC) environments to devices on the edge, i.e., physically near instruments of interest, has received tremendous interest in recent years. Such edge computing environments can operate on data in-situ, offering enticing benefits over data aggregation to HPC and CC facilities that include avoiding costs of transmission, increased data privacy, and real-time data analysis. Because of the inherent unreliability of edge computing environments, new fault tolerant approaches must be developed before the benefits of edge computing can be realized. Motivated by algorithm-based fault tolerance, a variant of the asynchronous Jacobi (ASJ) method is developed that achieves resilience to data corruption by rejecting solution approximations from neighbor devices according to a bound derived from convergence theory. Numerical results on a two-dimensional Poisson problem show the new rejection criterion, along with a novel approximation to the shortest path length on which the criterion depends, restores convergence for the ASJ variant in the presence of certain types data corruption. Numerical results are obtained for when the singular values in the analytic bound are approximated. A linear system with a more dense sparsity pattern is also explored. All results indicate that successful resilience to data corruption depends on whether the bound tightens fast enough to reject corrupted data before the iteration evolution deviates significantly from that predicted by the convergence theory defining the bound. This observation generalizes to future work on algorithm-based fault tolerance for other asynchronous algorithms, including upcoming approaches that leverage Krylov subspaces., Comment: Submitted to SIAM SISC September 29, 2023 Revision submitted to SIAM SISC April 7, 2024

Maintaining conservation laws in the fully discrete setting is critical for accurate long-time behavior of numerical simulations and requires accounting for discrete conservation properties in both space and time. This paper derives arbitrary order finite element exterior calculus spatial discretizations for the two-dimensional (2D) Navier-Stokes and drift-reduced magnetohydrodynamic equations that conserve both energy and enstrophy to machine precision when coupled with generally symplectic time-integration methods. Both continuous and discontinuous-Galerkin (DG) weak formulations can ensure conservation, but only generally symplectic time integration methods, such as the implicit midpoint method, permit exact conservation in time. Moreover, the symplectic implicit midpoint method yields an order of magnitude speedup over explicit schemes. The methods are implemented using the MFEM library and the solutions are verified for an extensive suite of 2D neutral fluid turbulence test problems. Numerical solutions are verified via comparison to a semi-analytic linear eigensolver as well as to the finite difference Global Drift Ballooning (GDB) code. However, it is found that turbulent simulations that conserve both energy and enstrophy tend to have too much power at high wavenumber and that this part of the spectrum should be controlled by reintroducing artificial dissipation. The DG formulation allows upwinding of the advection operator which dissipates enstrophy while still maintaining conservation of energy. Coupling upwinded DG with implicit symplectic integration appears to offer the best compromise of allowing mid-range wavenumbers to reach the appropriate amplitude while still controlling the high-wavenumber part of the spectrum.

As photovoltaic (PV) penetration continues to rise and smart inverter functionality continues to expand, smart inverters and other distributed energy resources (DERs) will play increasingly important roles in distribution system power management and security. In this paper, it is demonstrated that a constellation of smart inverters in a simulated distribution circuit can enable precise voltage predictions using an asynchronous and decentralized prediction algorithm. Using simulated data and a constellation of 15 inverters in a ring communication topology, the COLA algorithm is shown to accomplish the learning task required for voltage magnitude prediction with far less communication overhead than fully connected P2P learning protocols. Additionally, a dynamic stopping criterion is proposed that does not require a regularizer like the original COLA stopping criterion., Comment: This work has been submitted to the IEEE for possible publication

We analyze the stability and accuracy (up to third order) of a new family of implicit-explicit Runge-Kutta (IMEX RK) methods. This analysis expedites development of methods with various balances in the number of explicit stages and implicit solves. We emphasize deriving methods with large stability regions for horizontally explicit vertically implicit (HEVI) partitionings of nonhydrostatic atmosphere models. The IMKG2 and IMKG3 families of IMEX RK methods are formulated within this framework. The HOMME-NH model with a HEVI partitioning is used for testing the accuracy and stability of various IMKG2-3 methods. The efficiency of several IMKG2-3 methods is demonstrated in HOMME-NH and compared to other IMEX RK methods in the literature.

The nonhydrostatic High Order Method Modeling Environment (HOMME-NH) atmospheric dynamical core supports acoustic waves that propagate significantly faster than the advective wind speed, thus greatly limiting the timestep size that can be used with standard explicit time-integration methods. Resolving acoustic waves is unnecessary for accurate climate and weather prediction. This numerical stiffness is addressed herein by considering implicit-explicit additive Runge-Kutta (ARK IMEX) methods that can treat the acoustic waves in a stable manner without requiring implicit treatment of non-stiff modes. Various ARK IMEX methods are evaluated for their efficiency in producing accurate solutions, ability to take large timestep sizes, and sensitivity to grid cell length ratio. Both the Gravity Wave test and Baroclinic Instability test from the 2012 Dynamical Core Model Intercomparison Project (DCMIP) are used to recommend 5 of the 27 ARK IMEX methods for use in HOMME-NH., Comment: 22 pages, 3 Figures, 7 Tables

A high-resolution finite volume method approach to incorporating time-dependent slip across rectangular subfaults when modeling general fault geometry is presented. The fault slip is induced by a modification of the Riemann problem to the linear elasticity equations across cell interfaces aligned with the subfaults. This is illustrated in the context of the high-resolution wave-propagation algorithms that are implemented in the open source Clawpack software (www.clawpack.org), but this approach could be easily incorporated into other Riemann solver based numerical methods. Surface deformation results are obtained in both two and three dimensions and compared to those given by the steady-state, homogeneous half-space Okada solution.

The Closest Point method, initially developed by Ruuth and Merriman, allows for the numerical solution of surface partial differential equations without the need for a parameterization of the surface itself. Surface quantities are embedded into the surrounding domain by assigning each value at a given spatial location to the corresponding value at the closest point on the surface. This embedding allows for surface derivatives to be replaced by their Cartesian counterparts (e.g. $\nabla_s = \nabla$). This equivalence is only valid on the surface, and thus, interpolation is used to enforce what is known as the side condition away from the surface. To improve upon the method, this work derives an operator embedding that incorporates curvature information, making it valid in a neighborhood of the surface. With this, direct enforcement of the side condition is no longer needed. Comparisons in $\mathbb{R}^2$ and $\mathbb{R}^3$ show that the resulting Curvature-Augmented Closest Point method has better accuracy and requires less memory, through increased matrix sparsity, than the Closest Point method, while maintaining similar matrix condition numbers. To demonstrate the utility of the method in a physical application, simulations of inextensible, bi-lipid vesicles evolving toward equilibrium shapes are also included.