Search

We describe the ARKODE library of one-step time integration methods for ordinary differential equation (ODE) initial-value problems (IVPs). In addition to providing standard explicit and diagonally implicit Runge–Kutta methods, ARKODE also supports one-step methods designed to treat additive splittings of the IVP, including implicit-explicit (ImEx) additive Runge–Kutta methods and multirate infinitesimal (MRI) methods. We present the role of ARKODE within the SUNDIALS suite of time integration and nonlinear solver libraries, the core ARKODE infrastructure for utilities common to large classes of one-step methods, as well as its use of “time stepper” modules enabling easy incorporation of novel algorithms into the library. Numerical results show example problems of increasing complexity, highlighting the algorithmic flexibility afforded through this infrastructure, and include a larger multiphysics application leveraging multiple algorithmic features from ARKODE and SUNDIALS.

Recent studies have shown that reducing the precision of floating‐point calculations in an atmospheric model can improve the model's computational performance without affecting model fidelity, but code changes are needed to accommodate lower precision or to prevent undue round‐off error. For complex and computationally expensive systems like the Energy Exascale Earth System Model, a method is needed to objectively and efficiently assess the quality of the lower‐precision simulations and to quickly identify problematic code pieces. This paper demonstrates that the accuracy of the computed solution can be evaluated through a simple and quantitative error metric based on time step convergence. The proposed test method can unambiguously detect the accumulation of precision error and objectively assess its impact on solution accuracy. Furthermore, since fast physical processes are known to affect key features of the multiyear mean climate in atmospheric models, we show that the convergence test applied to short simulations can provide useful information about the impact of reduced precision on a model's long‐term behavior. In contrast, the traditional method of climate model evaluation requires multiple years of simulations and involves inspecting many physical quantities for statistical significance in the presence of natural variability. The simplicity, computational efficiency, and objectivity of the proposed method makes it very attractive for high‐resolution model development. The method is also expected to be applicable to other models that numerically solve time evolution equations.

Convergence testing is a common practice in the development of dynamical cores of atmospheric models but is not as often exercised for the parameterization of subgrid physics. An earlier study revealed that the stratiform cloud parameterizations in several predecessors of the Energy Exascale Earth System Model (E3SM) showed strong time step sensitivity and slower‐than‐expected convergence when the model's time step was systematically refined. In this work, a simplified atmosphere model is configured that consists of the spectral‐element dynamical core of the E3SM atmosphere model coupled with a large‐scale condensation parameterization based on commonly used assumptions. This simplified model also resembles E3SM and its predecessors in the numerical implementation of process coupling and shows poor time step convergence in short ensemble tests. We present a formal error analysis to reveal the expected time step convergence rate and the conditions for obtaining such convergence. Numerical experiments are conducted to investigate the root causes of convergence problems. We show that revisions in the process coupling and closure assumption help to improve convergence in short simulations using the simplified model; the same revisions applied to a full atmosphere model lead to significant changes in the simulated long‐term climate. This work demonstrates that causes of convergence issues in atmospheric simulations can be understood by combining analyses from physical and mathematical perspectives. Addressing convergence issues can help to obtain a discrete model that is more consistent with the intended representation of the physical phenomena.

Global atmospheric models seek to capture physical phenomena across a wide range of time and length scales. For this to be a feasible task, the physical processes with time or length scales below that of a computational time step or grid cell size are simplified as one or more parameterizations. Inadvertent oversimplification can violate constraints or destroy relationships in the original physical system and consequently lead to unexpected and physically invalid behavior. An example of such a problem has been investigated in the work of Wan et al. (2020, https://doi.org/10.1029/2019MS001982). This work addresses the issues at a more fundamental level by revisiting the parameterization derivation. A derivation of an unaveraged condensation rate in the unaveraged equations, sometimes referred to as subgrid equations, provides a clear description and more accurate quantification of the condensation/evaporation processes associated with cloud growth/decay, while avoiding simplifications used in earlier studies. A subgrid reconstruction (SGR) methodology is used to connect the unaveraged condensation rate with the grid cell averaged equations solved by the global model. Analyses of the SGR method and the numerical results provide insights into root causes of inconsistent discrete formulations and nonphysical behavior. It is also shown that the SGR methodology provides a flexible framework for addressing such inconsistencies. This work serves as a demonstration that when nonphysical behavior in a parameterization of subgrid variability is avoided through rigorous mathematical derivation, the resulting formulation can exhibit both better numerical convergence properties and significant impact on long‐term climate.

Transient simulation of power systems can be challenging. As models become more complex, applications require acceleration or faster than real time operations to assist system operators, and new methods of solving the problem are required. One potential approach includes time-parallel multigrid reduction methods that distribute the simulation time domain onto multiple compute units and solve the system equations with a parallel, iterative nonlinear multigrid scheme. Achieving significant speedup via these methods for real world power system problems requires examination of how the different features and characteristics of power system simulations interact with the parallel-in-time calculations and how they can be handled by such tools. One such characteristic is the use of unscheduled events that change the system at solution-dependent times. Examples of such events include equipment limits, governor deadbands, voltage regulator adjustments, or protection system actions. In sequential simulations, such events are detected by using root-finding methods or directly in control system models. In time-parallel methods these roots must be handled in parallel. In this paper we propose a method for addressing unscheduled events in time-parallel calculations for power system transient simulation.

The efficient simulation of non-hydrostatic atmospheric dynamics requires time integration methods capable of overcoming the explicit stability constraints on time step size arising from acoustic waves. In this work, we investigate various implicit–explicit (IMEX) additive Runge–Kutta (ARK) methods for evolving acoustic waves implicitly to enable larger time step sizes in a global non-hydrostatic atmospheric model. The IMEX formulations considered include horizontally explicit – vertically implicit (HEVI) approaches as well as splittings that treat some horizontal dynamics implicitly. In each case, the impact of solving nonlinear systems in each implicit ARK stage in a linearly implicit fashion is also explored.The accuracy and efficiency of the IMEX splittings, ARK methods, and solver options are evaluated on a gravity wave and baroclinic wave test case. HEVI splittings that treat some vertical dynamics explicitly do not show a benefit in solution quality or run time over the most implicit HEVI formulation. While splittings that implicitly evolve some horizontal dynamics increase the maximum stable step size of a method, the gains are insufficient to overcome the additional cost of solving a globally coupled system. Solving implicit stage systems in a linearly implicit manner limits the solver cost but this is offset by a reduction in step size to achieve the desired accuracy for some methods. Overall, the third-order ARS343 and ARK324 methods performed the best, followed by the second-order ARS232 and ARK232 methods.

Explicit Runge–Kutta methods and implicit multistep methods utilizing a Newton–Krylov nonlinear solver are evaluated for a range of configurations of the shallow-water dynamical core of the spectral element community atmosphere model to evaluate their computational performance. These configurations are designed to explore the attributes of each method under different but relevant model usage scenarios including varied spectral order within an element, static regional refinement, and scaling to the largest problem sizes. This analysis is performed within the shallow-water dynamical core option of a full climate model code base to enable a wealth of simulations for study, with the aim of informing solver development within the more complete hydrostatic dynamical core used for climate research. The limitations and benefits to using explicit versus implicit methods, with different parameters and settings, are discussed in light of the trade-offs with Message Passing Interface (MPI) communication and memory and their inherent efficiency bottlenecks. Given the performance behavior across the configurations analyzed here, the recommendation for future work using the implicit solvers is conditional based on scale separation and the stiffness of the problem. For the regionally refined configurations, the implicit method has about the same efficiency as the explicit method, without considering efficiency gains from a preconditioner. The potential for improvement using a preconditioner is greatest for higher spectral order configurations, where more work is shifted to the linear solver. Initial simulations with OpenACC directives to utilize a Graphics Processing Unit (GPU) when performing function evaluations show improvements locally, and that overall gains are possible with adjustments to data exchanges.

As part of the Exascale Computing Project (ECP), a recent focus of development efforts for the SUite of Nonlinear and DIfferential/ALgebraic equation Solvers (SUNDIALS) has been to enable GPU-accelerated time integration in scientific applications at extreme scales. This effort has resulted in several new GPU-enabled implementations of core SUNDIALS data structures, support for programming paradigms which are aware of the heterogeneous architectures, and the introduction of utilities to provide new points of flexibility. In this paper, we discuss our considerations, both internal and external, when designing these new features and present the features themselves. We also present performance results for several of the features on the Summit supercomputer and early access hardware for the Frontier supercomputer, which demonstrate negligible performance overhead resulting from the additional infrastructure and significant speedups when using both NVIDIA and AMD GPUs.

Productivity from day one on supercomputers that leverage new technologies requires significant preparation. An institution that procures a novel system architecture often lacks sufficient institutional knowledge and skills to prepare for it. Thus, the "Center of Excellence" (CoE) concept has emerged to prepare for systems such as Summit and Sierra, currently the top two systems in the Top 500. This paper documents CoE experiences that prepared a workload of diverse applications and math libraries for a heterogeneous system. We describe our approach to this preparation, including our management and execution strategies, and detail our experiences with and reasons for using different programming approaches. Our early science and performance results show that the project enabled significant early seismic science with up to a l4X throughput increase over Cori. In addition to our successes, we discuss our challenges and failures so others may benefit from our experience.

In this report we document performance test results on a SUNDIALS-based multiphysics demonstration application. We aim to assess the large-scale parallel performance of new capabilities that have been added to the SUNDIALS suite of time integrators and nonlinear solvers in recent years under funding from both the Exascale Computing Project (ECP) and the Scientific Discovery through Advanced Scientific (SciDAC) program, specifically: (a) SUNDIALS' new MPIManyVector module, that allows extreme flexibility in how a solution "vector" is staged on computational resources, (b) ARKode's new multirate integration module, MRIStep, allowing high-order accurate calculations that subcycle "fast" processes within "slow" ones, (c) SUNDIALS' new flexible linear solver interfaces, that allow streamlined specification of problem-specific linear solvers, and (d) SUNDIALS' new N_Vector additions of "fused" vector operations (to increase arithmetic intensity) and separation of reduction operations into "local" and "global" versions (to reduce latency by combining multiple reductions into a single MPI_Allreduce call). We anticipate that subsequent reports will extend this work to investigate a variety of other new features, including SUNDIALS' generic SUNNonlinearSolver interface and accelerator-enabled N_Vector modules, and upcoming MRIStep extensions to support custom "fast" integrators (that leverage problem structure) and IMEX integration of the "slow" time scale (to add diffusion).

Surface and subsurface flow constitute a naturally linked hydrologic continuum that has not traditionally been simulated in an integrated fashion. Recognizing the interactions between these systems has encouraged the development of integrated hydrologic models (IHMs) capable of treating surface and subsurface systems as a single integrated resource. IHMs is dynamically evolving with improvement in technology and the extent of their current capabilities are often only known to the developers and not general users. This article provides an overview of the core functionality, capability, applications, and ongoing development of one open-source IHM, ParFlow. ParFlow is a parallel, integrated, hydrologic model that simulates surface and subsurface flows. ParFlow solves Richards’ equation for three-dimensional variably saturated groundwater flow and the two-dimensional kinematic wave approximation of the shallow water equations for overland flow. The model employs a conservative centered finite difference scheme and a conservative finite volume method for subsurface flow and transport, respectively. ParFlow uses multigrid preconditioned Krylov and Newton-Krylov methods to solve the linear and nonlinear systems within each time step of the flow simulations. The code has demonstrated very efficient parallel solution capabilities. ParFlow has been coupled to geochemical reaction, land surface (e.g. Common Land Model), and atmospheric models to study the interactions among the subsurface, land surface, and the atmosphere systems across different spatial scales. This overview focuses on the current capabilities of the code, the core simulation engine, and the primary couplings of the subsurface model to other codes, taking a high-level perspective.

This work considers multirate generalized-structure additively partitioned Runge-Kutta (MrGARK) methods for solving stiff systems of ordinary differential equations (ODEs) with multiple time scales. These methods treat different partitions of the system with different timesteps for a more targeted and efficient solution compared to monolithic single rate approaches. With implicit methods used across all partitions, methods must find a balance between stability and the cost of solving nonlinear equations for the stages. In order to characterize this important trade-off, we explore multirate coupling strategies, problems for assessing linear stability, and techniques to efficiently implement Newton iterations for stage equations. Unlike much of the existing multirate stability analysis which is limited in scope to particular methods, we present general statements on stability and describe fundamental limitations for certain types of multirate schemes. New implicit multirate methods up to fourth order are derived, and their accuracy and efficiency properties are verified with numerical tests.

The efficient simulation of non-hydrostatic atmospheric dynamics requires time integration methods capable of overcoming the explicit stability constraints on time step size arising from acoustic waves. In this work we investigate various implicit-explicit (IMEX) additive Runge-Kutta (ARK) methods for evolving acoustic waves implicitly to enable larger time step sizes in a global non-hydrostatic atmospheric model. The IMEX formulations considered include horizontally implicit-vertically implicit (HEVI) approaches as well as splittings that treat some horizontal dynamics implicitly. In each case the impact of solving nonlinear systems in each implicit ARK stage in a linearly implicit fashion is also explored. The accuracy and efficiency of the IMEX splittings, ARK methods, and solver options are evaluated on a gravity wave and baroclinic wave test case. HEVI splittings that treat some vertical dynamics explicitly do not show a benefit in solution quality or run time over the most implicit HEVI formulation. While splittings that implicitly evolve some horizontal dynamics increase the maximum stable step size of a method, the gains are insufficient to overcome the additional cost of solving a globally coupled system. Solving implicit stage systems in a linearly implicit manner limits the solver cost but this is offset by a reduction in step size to achieve the desired accuracy. Overall, the third order ARS343 and ARK324 methods performed the best, followed by the second order ARS232 and ARK232 methods.

Due to complex dynamics inherent in the physical models, numerical formulation of subsurface and overland flow coupling can be challenging to solve. ParFlow is a subsurface flow code that utilizes a structured grid discretization in order to benefit from fast and efficient structured solvers. Implicit coupling between subsurface and overland flow modes in ParFlow is obtained by prescribing an overland boundary condition at the top surface of the computational domain. This form of implicit coupling leads to the activation and deactivation of the overland boundary condition during simulations where ponding or drying events occur. This results in a discontinuity in the discrete system that can be challenging to resolve. Furthermore, the coupling relies on unstructured connectivities between the subsurface and surface components of the discrete system, which makes it challenging to use structured solvers to effectively capture the dynamics of the coupled flow. We present a formulation of the discretized algebraic system that enables the use of an analytic form of the Jacobian for the Newton–Krylov solver, while preserving the structured properties of the discretization. An effective multigrid preconditioner is extracted from the analytic Jacobian and used to precondition the Jacobian linear system solver. We compare the performance of the new solver against one that uses a finite difference approximation to the Jacobian within the Newton–Krylov approach, previously used in the literature. Numerical results explores the effectiveness of using the analytic Jacobian for the Newton–Krylov solver, and highlights the performance of the new preconditioner and its cost. The results indicate that the new solver is robust and generally outperforms the solver that is based on the finite difference approximation to the Jacobian, for problems where the overland boundary condition is activated and deactivated during the simulation. A parallel weak scaling study highlights the efficiency of the new solver.

Performing accurate and efficient numerical simulation of global atmospheric climate models is challenging due to the disparate length and time scales over which physical processes interact. Implicit solvers enable the physical system to be integrated with a time step commensurate with processes being studied. The dominant cost of an implicit time step is the ancillary linear system solves, so we have developed a preconditioner aimed at improving the efficiency of these linear system solves. Our preconditioner is based on an approximate block factorization of the linearized shallow-water equations and has been implemented within the spectral element dynamical core within the Community Atmospheric Model (CAM-SE). In this paper, we discuss the development and scalability of the preconditioner for a suite of test cases with the implicit shallow-water solver within CAM-SE.

Multiphysics simulations frequently are composed from highly-optimized solvers for physical subprocesses through the use of temporal operator splitting. The subphysics are evolved sequentially, passing information between solver components as needed. It is often useful yet difficult to determine how the simulation error depends on the temporal splitting method, the time step size and the discretization parameters in the solver components. This paper proposes a framework to decompose the total error in a quantity of interest for an advection–diffusion simulation into two primary contributions: that due to the operator splitting, as implied by the structure of the code, and that due to the discretization errors from the component solvers. The method is applied to the advection–diffusion equation with boundaries. The advection and diffusion operators require separate boundary conditions upon splitting; the specification of these impacts the splitting error contribution. Computational examples demonstrate that the proposed method successfully identifies the splitting error contribution, including that induced by a suboptimal imposition of boundary conditions, and the decrease in the splitting error contribution in moving to a higher-order splitting method. The discretization error contribution is decomposed further into contributions attributed to the separate advection and diffusion solvers.

This paper presents a multigrid reduction in time (MGRIT) algorithm for achieving time parallelism using multistep backward difference formula (BDF) methods on variably-spaced temporal grids. This MGRIT approach transforms the linear multistep methods into single step methods applied to groups of time steps. Stability considerations are addressed through lowering of the order on coarse grids. The methods are presented for fixed and variable time step formulations. Numerical results show moderate speedups for the heat equation as well as two IEEE power grid test problems characterized by nonlinear differential algebraic equations (DAE) of index 1. The performance of MGRIT with BDF is compared to MGRIT with Runge-Kutta for both the fixed and variable step methods of the same order.

We consider multiphysics applications from algorithmic and architectural perspectives, where “algorithmic” includes both mathematical analysis and computational complexity, and “architectural” includes both software and hardware environments. Many diverse multiphysics applications can be reduced, en route to their computational simulation, to a common algebraic coupling paradigm. Mathematical analysis of multiphysics coupling in this form is not always practical for realistic applications, but model problems representative of applications discussed herein can provide insight. A variety of software frameworks for multiphysics applications have been constructed and refined within disciplinary communities and executed on leading-edge computer systems. We examine several of these, expose some commonalities among them, and attempt to extrapolate best practices to future systems. From our study, we summarize challenges and forecast opportunities.

This paper is concerned with the computation of numerical discretization error for uncertainty quantification. An a posteriori error formula is described for a functional measurement of the solution to a scalar advection equation that is estimated by finite volume approximations. An exact error formula and computable error estimate are derived based on an abstractly defined approximation of the adjoint solution. The adjoint problem is divorced from the finite volume method used to approximate the forward solution variables and may be approximated using a low-order finite volume method. The accuracy of the computable error estimate provably satisfies an a priori error bound for sufficiently smooth solutions of the forward and adjoint problems. Computational examples are provided that show support of the theory for smooth solutions. The application to problems with discontinuities is also investigated computationally.

This paper presents a fully multilevel approach to parallel in time solution of transient power system simulations. The method employs a multigrid reduction algorithm in time parallelized using the MPI distributed memory programming model. The method is demonstrated on a simple Single Machine Infinite Bus differential-algebraic equation model problem, for which speedup is obtained for as few as 8 processing cores on a problem with 10,000 time steps. Speedup of a factor of 13 is observed for a 100,000 step version of this simple problem. Based on these results, we expect significantly better speedup on larger problems where more work is available to each processor allowing greater amortization of the parallel communication.

The estimation of discretization error in numerical simulations is a key component in the development of uncertainty quantification. In particular, there exists a need for reliable, robust estimators for finite volume and finite difference discretizations of hyperbolic partial differential equations. The approach espoused here, often called the error transport approach in the literature, is to solve an auxiliary error equation concurrently with the primal governing equation to obtain a point-wise (cell-wise) estimate of the discretization error. Nonlinear, time-dependent problems are considered. In contrast to previous work, fully nonlinear error equations are advanced, and potential benefits are identified. A systematic approach to approximate the local residual for both method-of-lines and space–time discretizations is developed. Behavior of the error estimates on problems that include weak solutions demonstrates the positive properties of nonlinear error transport.

Complete models of the hydrologic cycle have gained recent attention as research has shown interdependence between the coupled land and energy balance of the subsurface, land surface, and lower atmosphere. PF.WRF is a new model that is a combination of the Weather Research and Forecasting (WRF) atmospheric model and a parallel hydrology model (ParFlow) that fully integrates three-dimensional, variably saturated subsurface flow with overland flow. These models are coupled in an explicit, operator-splitting manner via the Noah land surface model (LSM). Here, the coupled model formulation and equations are presented and a balance of water between the subsurface, land surface, and atmosphere is verified. The improvement in important physical processes afforded by the coupled model using a number of semi-idealized simulations over the Little Washita watershed in the southern Great Plains is demonstrated. These simulations are initialized with a set of offline spinups to achieve a balanced state of initial conditions. To quantify the significance of subsurface physics, compared with other physical processes calculated in WRF, these simulations are carried out with two different surface spinups and three different microphysics parameterizations in WRF. These simulations illustrate enhancements to coupled model physics for two applications: water resources and wind-energy forecasting. For the water resources example, it is demonstrated how PF.WRF simulates explicit rainfall and water storage within the basin and runoff. Then the hydrographs predicted by different microphysics schemes within WRF are compared. Because soil moisture is expected to impact boundary layer winds, the applicability of the model to wind-energy applications is demonstrated by using PF.WRF and WRF simulations to provide estimates of wind and wind shear that are useful indicators of wind-power output.

We introduce an operator-based scheme for preconditioning stiff components encountered in implicit methods for hyperbolic systems of PDEs posed on regular grids. The method is based on a directional splitting of the implicit operator, followed by a characteristic decomposition of the resulting directional parts. This approach allows for the solution of any number of characteristic components, from the entire system to only the fastest, stiffness-inducing waves. We apply the preconditioning method to stiff hyperbolic systems arising in magnetohydrodynamics and gas dynamics. We then present numerical results showing that this preconditioning scheme works well on problems where the underlying stiffness results from the interaction of fast transient waves with slowly-evolving dynamics, scales well to large problem sizes and numbers of processors, and allows for additional customization based on the specific problems under study.

A Newton-Krylov method is an implementation of Newton's method in which a Krylov subspace method is used to solve approximately the linear systems that characterize steps of Newton's method. Newton-Krylov methods are often implemented in “matrix-free” form, in which the Jacobian-vector products required by the Krylov solver are approximated by finite differences. Here we consider using approximate function values in these finite differences. We first formulate a finite-difference Arnoldi process that uses approximate function values. We then outline a Newton-Krylov method that uses an implementation of the GMRES or Arnoldi method based on this process, and we develop a local convergence analysis for it, giving sufficient conditions on the approximate function values for desirable local convergence properties to hold. We conclude with numerical experiments involving particular function-value approximations suitable for nonlinear diffusion problems. For this case, conditions are given for meeting the convergence assumptions for both lagging and linearizing the nonlinearity in the function evaluation.

Recent work suggests that Anderson acceleration can be used as an accelerator to the fixed-point iterative method. To improve the viability of the algorithm, we seek to improve its computational efficiency on parallel machines. The primary kernel of the method is a least-squares minimization within the main loop. We consider two approaches to reduce its cost. The first is to use a communication-avoiding QR factorization, and the second is to employ a GMRES-like restarting procedure. On problems using 1,000 processors or less, we find the amount of communication too low to justify communication avoidance. The restarting procedure also proves not to be better than current approaches unless the cost of the function evaluation is very small. In order to begin taking advantage of current trends in machine architecture, we also studied a first-attempt single-node GPU implementation of Anderson acceleration. Performance results show that for sufficiently large problems a GPU implementation can provide a significant performance increase over CPU versions due to the GPU’s higher memory bandwidth.

SUNDIALS is a suite of advanced computational codes for solving large-scale problems that can be modeled as a system of nonlinear algebraic equations, or as initial-value problems in ordinary differential or differential-algebraic equations. The basic versions of these codes are called KINSOL, CVODE, and IDA, respectively. The codes are written in ANSI standard C and are suitable for either serial or parallel machine environments. Common and notable features of these codes include inexact Newton-Krylov methods for solving large-scale nonlinear systems; linear multistep methods for time-dependent problems; a highly modular structure to allow incorporation of different preconditioning and/or linear solver methods; and clear interfaces allowing for users to provide their own data structures underneath the solvers. We describe the current capabilities of the codes, along with some of the algorithms and heuristics used to achieve efficiency and robustness. We also describe how the codes stem from previous and widely used Fortran 77 solvers, and how the codes have been augmented with forward and adjoint methods for carrying out first-order sensitivity analysis with respect to model parameters or initial conditions.

© The Authors. Published by Elsevier B.V. We present a survey of physical and computational techniques that have the potential to contribute to the next generation of high-fidelity, multi-scale climate simulations. Examples of the climate science problems that can be investigated with more depth with these computational improvements include the capture of remote forcings of localized hydrological extreme events, an accurate representation of cloud features over a range of spatial and temporal scales, and parallel, large ensembles of simulations to more effectively explore model sensitivities and uncertainties. Numerical techniques, such as adaptive mesh refinement, implicit time integration, and separate treatment of fast physical time scales are enabling improved accuracy and fidelity in simulation of dynamics and allowing more complete representations of climate features at the global scale. At the same time, partnerships with computer science teams have focused on taking advantage of evolving computer architectures such as many-core processors and GPUs. As a result, approaches which were previously considered prohibitively costly have become both more efficient and scalable. In combination, progress in these three critical areas is poised to transform climate modeling in the coming decades. These topics have been presented within a workshop titled, "Numerical and Computational Developments to Advance Multiscale Earth System Models (MSESM'15)," as part of the International Conference on Computational Sciences, Reykjavik, Iceland, June 1-3, 2015.

We present an analysis of expanded mixed finite element methods applied to Richards' equation, a nonlinear parabolic partial differential equation modeling the flow of water into a variably saturated porous medium. We consider the full range of saturated to completely unsaturated media. In the case of the lowest order Raviart--Thomas spaces and the range of all possible saturations, we bound the H-1-norm of the error in capacity in terms of approximation error. This estimate uses a time-integrated scheme and the Kirchhoff transformation to handle a degeneracy in the case of completely unsaturated flow. Optimal convergence is then shown for a nonlinear form related to the error in the capacity for the case of saturated to partially saturated flow. Convergence rates depending on the Holder continuity of the capacity term are derived. Last, optimal convergence of pressures and fluxes is stated for the case of strictly partially saturated flow.

We present a two-level finite difference scheme for the approximation of nonlinear parabolic equations. Discrete inner products and the lowest-order Raviart--Thomas approximating space are used in the expanded mixed method in order to develop the finite difference scheme. Analysis of the scheme is given assuming an implicit time discretization. In this two-level scheme, the full nonlinear problem is solved on a "coarse" grid of size H. The nonlinearities are expanded about the coarse grid solution and an appropriate interpolation operator is used to provide values of the coarse grid solution on the fine grid in terms of superconvergent node points. The resulting linear but nonsymmetric system is solved on a "fine" grid of size h. Some a priori error estimates are derived which show that the discrete L\infty(L2) and L2(H1) errors are $O(h^2 + H^{4-d/2} + \Delta t)$, where $d \geq 1$ is the spatial dimension.

An a posteriori error formula is described when a statistical measurement of the solution to a hyperbolic conservation law in 1D is estimated by finite volume approximations. This is accomplished using adjoint error estimation. In contrast to previously studied methods, the adjoint problem is divorced from the finite volume method used to approximate the forward solution variables. An exact error formula and computable error estimate are derived based on an abstractly defined approximation of the adjoint solution. This framework allows the error to be computed to an arbitrary accuracy given a sufficiently well resolved approximation of the adjoint solution. The accuracy of the computable error estimate provably satisfies an a priori error bound for sufficiently smooth solutions of the forward and adjoint problems. The theory does not currently account for discontinuities. Computational examples are provided that show support of the theory for smooth solutions. The application to problems with discontinuities is also investigated computationally.

[1] We present the results of a unique, parallel scaling study using a 3-D variably saturated flow problem including land surface processes that ranges from a single processor to a maximum number of 16,384 processors. In the applied finite difference framework and for a fixed problem size per processor, this results in a maximum number of approximately 8 × 109 grid cells (unknowns). Detailed timing information shows that the applied simulation platform ParFlow exhibits excellent parallel efficiency. This study demonstrates that regional scale hydrologic simulations on the order of 103 km2 are feasible at hydrologic resolution (∼100–101 m laterally, 10−2–10−1 m vertically) with reasonable computation times, which has been previously assumed to be an intractable computational problem.

In this paper, the authors investigate performance of a fully implicit formulation and solution method of a diffusion-reaction system modeling radiation diffusion with material energy transfer and a fusion fuel source. In certain parameter regimes this system can lead to a rapid conversion of potential energy into material energy. Accuracy in time integration is essential for a good solution since a major fraction of the fuel can be depleted in a very short time. Such systems arise in a number of application areas including evolution of a star and inertial confinement fusion. Previous work has addressed implicit solution of radiation diffusion problems. Recently Shadid and coauthors have looked at implicit and semi-implicit solution of reaction-diffusion systems. In general they have found that fully implicit is the most accurate method for difficult coupled nonlinear equations. In previous work, they have demonstrated that a method of lines approach coupled with a BDF time integrator and a Newton-Krylov nonlinear solver could efficiently and accurately solve a large-scale, implicit radiation diffusion problem. In this paper, they extend that work to include an additional heating term in the material energy equation and an equation to model the evolution of the reactive fuel density. This systemmore » now consists of three coupled equations for radiation energy, material energy, and fuel density. The radiation energy equation includes diffusion and energy exchange with material energy. The material energy equation includes reaction heating and exchange with radiation energy, and the fuel density equation includes its depletion due to the fuel consumption.« less

This paper provides extensions of an element agglomeration AMG method to nonlinear elliptic problems discretized by the finite element method on general unstructured meshes. The method constructs coarse discretization spaces and corresponding coarse nonlinear operators as well as their Jacobians. We introduce both standard (fairly quasi-uniformly coarsened) and non-standard (coarsened away) coarse meshes and respective finite element spaces. We use both kind of spaces in FAS type coarse subspace correction (or Schwarz) algorithms. Their performance is illustrated on a number of model problems. The coarsened away spaces seem to perform better than the standard spaces for problems with nonlinearities in the principal part of the elliptic operator.

In dislocation dynamics simulations, strain hardening simulations require integrating stiff systems of ordinary differential equations in time with expensive force calculations, discontinuous topological events and rapidly changing problem size. Current solvers in use often result in small time steps and long simulation times. Faster solvers may help dislocation dynamics simulations accumulate plastic strains at strain rates comparable to experimental observations. This paper investigates the viability of high-order implicit time integrators and robust nonlinear solvers to reduce simulation run times while maintaining the accuracy of the computed solution. In particular, implicit Runge–Kutta time integrators are explored as a way of providing greater accuracy over a larger time step than is typically done with the standard second-order trapezoidal method. In addition, both accelerated fixed point and Newton's method are investigated to provide fast and effective solves for the nonlinear systems that must be resolved within each time step. Results show that integrators of third order are the most effective, while accelerated fixed point and Newton's method both improve solver performance over the standard fixed point method used for the solution of the nonlinear systems.

In this paper, we present results demonstrating the effectiveness of a sensitivity analysis approach to uncertainty quantification of a variably saturated flow model. The basis for our method is a software system which simultaneously solves for solutions of large-scale nonlinear systems of equations and the sensitivity of the solutions to selected parameters. We present test cases showing the effects on the relative uncertainty of pressure due to heterogeneity in the absolute permeability and to differences in parameterizing the Van Genuchten curve soil parameters, α and n .

This paper describes an implicit approach and nonlinear solver for solution of radiation-hydrodynamic problems in the context of supernovae and proto-neutron star cooling. The robust approach applies Newton-Krylov methods and overcomes the difficulties of discontinuous limiters in the discretized equations and scaling of the equations over wide ranges of physical behavior. We discuss these difficulties, our approach for overcoming them, and numerical results demonstrating accuracy and efficiency of the method.

Efficient solutions of global climate models require effectively handling disparate length and time scales. Implicit solution approaches allow time integration of the physical system with a step size governed by accuracy of the processes of interest rather than by stability of the fastest time scales present. Implicit approaches, however, require the solution of nonlinear systems within each time step. Usually, a Newton's method is applied to solve these systems. Each iteration of the Newton's method, in turn, requires the solution of a linear model of the nonlinear system. This model employs the Jacobian of the problem-defining nonlinear residual, but this Jacobian can be costly to form. If a Krylov linear solver is used for the solution of the linear system, the action of the Jacobian matrix on a given vector is required. In the case of spectral element methods, the Jacobian is not calculated but only implemented through matrix- vector products. The matrix-vector multiply can also be approximated by a finite difference approximation which may introduce inaccuracy in the overall nonlinear solver. In this paper, we review the advantages and disadvantages of finite difference approximations of these matrix- vector products for climate dynamics within the spectral element shallow water dynamical core of the Community Atmosphere Model (CAM).