Search

Discretized numerical models of the atmosphere are usually intended to faithfully represent an underlying set of continuous equations, but this necessary condition is violated sometimes by subtle pathologies that have crept into the discretized equations. Such pathologies can introduce undesirable artifacts, such as sawtooth noise, into the model solutions. The presence of these pathologies can be detected by numerical convergence testing. This study employs convergence testing to verify the discretization of the Cloud Layers Unified By Binormals (CLUBB) model of clouds and turbulence. That convergence testing identifies two aspects of CLUBB’s equation set that contribute to undesirable noise in the solutions. First, numerical limiters (i.e. clipping) used by CLUBB introduce discontinuities or slope discontinuities in model fields. Second, nonlinear numerical diffusion employed for improving numerical stability can introduce unintended small-scale features into the solution of the model equations. Smoothing the limiters and using linear diffusion (low-order hyperdiffusion) reduces the noise and restores the expected first-order convergence in CLUBB’s solutions. These model reformulations enhance our confidence in the trustworthiness of solutions from CLUBB by eliminating the unphysical oscillations in high-resolution simulations. The improvements in the results at coarser, near-operational grid spacing and timestep are also seen in cumulus cloud and dry turbulence tests. In addition, convergence testing is proven to be a valuable tool for detecting pathologies, including unintended discontinuities and grid dependence, in the model equation set.

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.

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., 22 pages, 3 Figures, 7 Tables

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.

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.

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.

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.

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.

Many multiphysics applications arise in the world of mathematical modeling and simulation. Much of the time in scientific computation these multiphysics applications are solved by decoupling the physics, giving no heed to how this affects the numerical results. However, a fully coupled approach is often not computationally cost effective. Consequently, having a metric for determining the strength of coupling could give insight into whether a simulation should be decoupled in the computation. If the fully coupled approach is not available, then a metric that measures the strength of coupling dynamically in time could help determine when smaller time steps are required to better incorporate coupling into the split solution. In this paper, we report on an Institute for Mathematics and Its Applications student project where we explored metrics for dynamically measuring the strength of coupling between two physical components in a model multiphysics simulation. Four metrics were considered: two based on measured components of the Jacobian matrix, one on error estimates, and the last on timescales of the system components. The metrics are all developed based on the previous work found in the literature and tested on a diffusion–reaction problem.

In this paper, we investigate the effectiveness of the Anderson acceleration method applied to modified Picard iteration for nonlinear problems arising in variably saturated flow modeling. While many authors have studied the relative merits of Newton’s method and modified Picard iteration in this context, the combination of Anderson acceleration and modified Picard iteration has not been investigated for these problems until recently. Since modified Picard iteration can be slow to converge, we investigate the use of Anderson acceleration to provide faster convergence while maintaining the robustness and lower memory requirements of modified Picard iteration relative to Newton’s method. Results indicate that Anderson acceleration significantly improves not only convergence speed but also robustness of modified Picard iteration and can often provide faster solutions than Newton’s method without the need for derivative computations.

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.

This paper introduces a federated simulation toolkit (FSKIT) that couples continuous time and discrete event simula- tions (DES) to perform the co-simulation of electric power grids and communication networks. A High Performance Computing (HPC) oriented power system dynamic simulator, GridDyn, was used for the electric power grid simulation. GridDyn is coupled to the open-source network simulator, ns-3, through FSKIT. FSKIT provides time control for advancing the state of federated simulators, and facilitates communication among objects in the federate. A wide-area communication-based electric transmission protection scheme is simulated with FSKIT, using the IEEE 39- bus test system. A communication network for the 39-bus system is built in ns-3, and basic protection relay logic is added to the power system model in order to perform the co-simulation.

© 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 a nonlinearly implicit, conservative numerical method for integration of the single-fluid resistive MHD equations. The method uses a high-order spatial discretization that preserves the solenoidal property of the magnetic field. The fully coupled PDE system is solved implicitly in time, providing for increased interaction between physical processes as well as additional stability over explicit-time methods. A high-order adaptive time integration is employed, which in many cases enables time steps ranging from one to two orders of magnitude larger than those constrained by the explicit CFL condition. We apply the solution method to illustrative examples relevant to stiff magnetic fusion processes which challenge the efficiency of explicit methods. We provide computational evidence showing that for such problems the method is comparably accurate with explicit-time simulations, while providing a significant runtime improvement due to its increased temporal stability.

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.

We present a solution method for fully implicit radiation diffusion problems discretized on meshes having millions of spatial zones. This solution method makes use of high order in time integration techniques, inexact Newton-Krylov nonlinear solvers, and multigrid preconditioners. We explore the advantages and disadvantages of high order time integration methods for the fully implicit formulation on both two- and three-dimensional problems with tabulated opacities and highly nonlinear fusion source terms.

In this paper, we discuss the computation and use of solution sensitivities for analyzing radiation diffusion problems and the dependence of solutions on input parameters. The derivation of the sensitivity equations is given, along with a description of the method for solving them in tandem with the simulation. The parameter values express material opacity as a power-law of material temperature and density. The computed sensitivities reveal important qualitative details about the temperature coupling and diffusion processes. It is also shown that these sensitivities are valuable for ranking the parameters from most to least influential, designing improved experiments, and quantifying uncertainty in the simulation results. Lastly, the numerical examples show that these various types of sensitivity analysis are only moderately expensive to perform relative to solving the simulation by itself.

In this paper we revisit and prove optimal-order and mesh-independent convergence of an inexact Newton method where the linear Jacobian systems are solved with multigrid techniques. This convergence is shown using Banach spaces and the norm $\max\{\|\cdot\|_1,\; \|\cdot\|_{0,\infty}\}$, a stronger norm than is used in previous work. These results are valid for a class of second-order, semilinear, finite element, elliptic problems posed on quasi-uniform grids. Numerical results are given which validate the theory.

In this paper, parallelization and high performance computing are utilized to enable ultrafast transient stability analysis that can be used in a real-time environment to quickly perform “what-if” simulations involving system dynamics phenomena. EPRI's Extended Transient Midterm Simulation Program (ETMSP) is modified and enhanced for this work. The contingency analysis is scaled for large-scale contingency analysis using Message Passing Interface (MPI) based parallelization. Simulations of thousands of contingencies on a high performance computing machine are performed, and results show that parallelization over contingencies with MPI provides good scalability and computational gains. Different ways to reduce the Input/Output (I/O) bottleneck are explored, and findings indicate that architecting a machine with a larger local disk and maintaining a local file system significantly improve the scaling results. Thread-parallelization of the sparse linear solve is explored also through use of the SuperLU_MT library.

In this paper, we present a class of solvers developed for the parallel solution of Richards' equation, a model used in variably saturated flow simulations. These solvers take advantage of the fast, robust convergence of globalized Newton methods as well as the parallel scalability of multigrid preconditioners. We compare two multigrid methods. The methods differ primarily in their handling of discontinuous and anisotropic permeability fields, with one method invoking a simple pointwise smoothing technique and the other a more expensive plane smoother. Computational results are presented to show the effectiveness of the entire nonlinear solution procedure, to demonstrate the effect of discontinuities and anisotropies, and to explore parallel efficiencies.

In this paper, we present a comparison of four preconditioning strategies for Jacobian systems arising in the fully implicit solution of radiation diffusion coupled with material energy transfer. The four preconditioning methods are block Jacobi, Schur complement, and operator splitting approaches that split the preconditioner solve into two steps. One splitting method includes the coupling of the radiation and material fields that appears in the matrix diagonal in the first solve, and the other method puts this coupling into the second solve. All preconditioning approaches use multigrid methods to invert blocks of the matrix formed from the diffusion operator. The Schur complement approach is clearly seen to be the most effective for a large range of weightings between the diffusion and energy coupling terms. In addition, tabulated opacity studies were conducted where, again, the Schur preconditioner performed well. Last, a parallel scaling study was done showing algorithmic scalability of the Schur preconditioner.

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.

A new parallel solution technique is developed for the fully implicit three‐dimensional two‐phase flow model. An expandedcell‐centered finite difference scheme which allows for a full permeability tensor is employed for the spatial discretization, and backwardEuler is used for the time discretization. The discrete systems are solved using a novel inexact Newton method that reuses the Krylov information generated by the GMRES linear iterative solver. Fast nonlinear convergence can be achieved by composing inexact Newton steps with quasi‐Newton steps restricted to the underlying Krylov subspace. Furthermore, robustness and efficiency are achieved with a line‐search backtracking globalization strategy for the nonlinear systems and a preconditioner for each coupled linear system to be solved. This inexact Newton method also makes use of forcing terms suggested by Eisenstat and Walker which prevent oversolving of the Jacobian systems. The preconditioner is a new two‐stage method which involves a decoupling strategy plus the separate solutions of both nonwetting‐phase pressure and saturation equations. Numerical results show that these nonlinear and linear solvers are very effective.

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.