2018

U. Ruede, K. Willcox, L. C. McInnes, H. De Sterck, G. Biros, H. Bungartz, J. Corones, E. Cramer, J. Crowley, O. Ghattas, M. Gunzburger, M. Hanke, R. Harrison, M. Heroux, J. Hesthaven, P. Jimack, C. Johnson, K. E. Jordan, D. E. Keyes, R. Krause, V. Kumar, S. Mayer, J. Meza, K. M. Mrken, J. T. Oden, L. Petzold, P. Raghavan, S. M. Shontz, A. Trefethen, P. Turner, V. Voevodin, B. Wohlmuth,, C. S. Woodward.
**“Research and Education in Computational Science and Engineering,”** Subtitled **“Report from a workshop sponsored by the Society for Industrial and Applied Mathematics (SIAM) and the European Exascale Software Initiative (EESI-2), August 4-6, 2014, Breckenridge, Colorado,”** Vol. abs/1610.02608, 2018.

This report presents challenges, opportunities and directions for computational science and engineering (CSE) research and education for the next decade. Over the past two decades the field of CSE has penetrated both basic and applied research in academia, industry, and laboratories to advance discovery, optimize systems, support decision-makers, and educate the scientific and engineering workforce. Informed by centuries of theory and experiment, CSE performs computational experiments to answer questions that neither theory nor experiment alone is equipped to answer. CSE provides scientists and engineers with algorithmic inventions and software systems that transcend disciplines and scales. CSE brings the power of parallelism to bear on troves of data. Mathematics-based advanced computing has become a prevalent means of discovery and innovation in essentially all areas of science, engineering, technology, and society; and the CSE community is at the core of this transformation. However, a combination of disruptive developments—including the architectural complexity of extreme-scale computing, the data revolution and increased attention to data-driven discovery, and the specialization required to follow the applications to new frontiers—is redefining the scope and reach of the CSE endeavor. With these many current and expanding opportunities for the CSE field, there is a growing demand for CSE graduates and a need to expand CSE educational offerings. This need includes CSE programs at both the undergraduate and graduate levels, as well as continuing education and professional development programs, exploiting the synergy between computational science and data science. Yet, as institutions consider new and evolving educational programs, it is essential to consider the broader research challenges and opportunities that provide the context for CSE education and workforce development.

W.Thevathasan, B. Debu, T. Aziz, B. R. Bloem, C. Blahak, C. Butson, V. Czernecki, T. Foltynie, V. Fraix, D. Grabli, C. Joint, A. M. Lozano, M. S. Okun, J. Ostrem, N. Pavese, C. Schrader, C. H. Tai, J. K. Krauss, E. Moro.
**“Pedunculopontine nucleus deep brain stimulation in Parkinson's disease: A clinical review,”** In *Movement Disorders*, Vol. 33, No. 1, pp. 10--20. 2018.

ISSN: 1531-8257

DOI: 10.1002/mds.27098

Pedunculopontine nucleus region deep brain stimulation (DBS) is a promising but experimental therapy for axial motor deficits in Parkinson's disease (PD), particularly gait freezing and falls. Here, we summarise the clinical application and outcomes reported during the past 10 years. The published dataset is limited, comprising fewer than 100 cases. Furthermore, there is great variability in clinical methodology between and within surgical centers. The most common indication has been severe medication refractory gait freezing (often associated with postural instability). Some patients received lone pedunculopontine nucleus DBS (unilateral or bilateral) and some received costimulation of the subthalamic nucleus or internal pallidum. Both rostral and caudal pedunculopontine nucleus subregions have been targeted. However, the spread of stimulation and variance in targeting means that neighboring brain stem regions may be implicated in any response. Low stimulation frequencies are typically employed (20-80 Hertz). The fluctuating nature of gait freezing can confound programming and outcome assessments. Although firm conclusions cannot be drawn on therapeutic efficacy, the literature suggests that medication refractory gait freezing and falls can improve. The impact on postural instability is unclear. Most groups report a lack of benefit on gait or limb akinesia or dopaminergic medication requirements. The key question is whether pedunculopontine nucleus DBS can improve quality of life in PD. So far, the evidence supporting such an effect is minimal. Development of pedunculopontine nucleus DBS to become a reliable, established therapy would likely require a collaborative effort between experienced centres to clarify biomarkers predictive of response and the optimal clinical methodology.

2017

M. Berzins, D. A. Bonnell, Jr. Cizewski, K. M. Heeger, A.J.G. Hey, C. J. Keane, B. A. Ramsey, K. A. Remington, J.L. Rempe.
**“Department of Energy, Advanced Scientific Computing Advisory Committee (ASCAC), Subcommittee on LDRD Review Final Report,”** May, 2017.

M. Berzins.
**“Nonlinear Stability of the MPM Method,”** In *V International Conference on Particle-based Methods – Fundamentals and Applications. PARTICLES 2017*, Edited by P. Wriggers, M. Bischoff, E. O˜nate, D.R.J. Owen, & T. Zohdi, 2017.

The Material Point Method (MPM) has been very successful in providing solutions to many challenging problems involving large deformations. The nonlinear nature of MPM makes it necessary to use a full nonlinear stability analysis to determine a stable timestep. The stability analysis of Spigler and Vianello is adapted to MPM and used to derive a stable timestep bound for a model problem. This bound is contrasted against a traditional CFL bound.

A. Bhatele, J. Yeom, N. Jain, C. J. Kuhlman, Y. Livnat, K. R. Bisset, L. V. Kale, M. V. Marathe.
**“Massively Parallel Simulations of Spread of Infectious Diseases over Realistic Social Networks,”** In *2017 17th IEEE/ACM International Symposium on Cluster, Cloud and Grid Computing (CCGRID)*, May, 2017.

DOI: 10.1109/ccgrid.2017.141

Controlling the spread of infectious diseases in large populations is an important societal challenge. Mathematically, the problem is best captured as a certain class of reaction-diffusion processes (referred to as contagion processes) over appropriate synthesized interaction networks. Agent-based models have been successfully used in the recent past to study such contagion processes. We describe EpiSimdemics, a highly scalable, parallel code written in Charm++ that uses agent-based modeling to simulate disease spreads over large, realistic, co-evolving interaction networks. We present a new parallel implementation of EpiSimdemics that achieves unprecedented strong and weak scaling on different architectures — Blue Waters, Cori and Mira. EpiSimdemics achieves five times greater speedup than the second fastest parallel code in this field. This unprecedented scaling is an important step to support the long term vision of real-time epidemic science. Finally, we demonstrate the capabilities of EpiSimdemics by simulating the spread of influenza over a realistic synthetic social contact network spanning the continental United States (∼280 million nodes and 5.8 billion social contacts).

L. Bos, A. Narayan, N. Levenberg, F. Piazzon.
**“An Orthogonality Property of the Legendre Polynomials,”** In *Constructive Approximation*, Vol. 45, No. 1, pp. 65--81. Feb, 2017.

ISSN: 0176-4276, 1432-0940

DOI: 10.1007/s00365-015-9321-3

We give a remarkable additional othogonality property of the classical Legendre polynomials on the real interval [−1,1]: polynomials up to degree n from this family are mutually orthogonal under the arcsine measure weighted by the degree-n normalized Christoffel function

J. Cates, L. Nevell, S. I. Prajapati, L. D. Nelon, J. Y. Chang, M. E. Randolph, B. Wood, C. Keller, R. T. Whitaker.
**“Shape analysis of the basioccipital bone in Pax7-deficient mice,”** In *Scientific Reports*, Vol. 7, No. 1, Springer Nature, Dec, 2017.

DOI: 10.1038/s41598-017-18199-9

We compared the cranial base of newborn *Pax7*-deficient and wildtype mice using a computational shape modeling technology called particle-based modeling (PBM). We found systematic differences in the morphology of the basiooccipital bone, including a broadening of the basioccipital bone and an antero-inferior inflection of its posterior edge in the *Pax7*-deficient mice. We show that the *Pax7* cell lineage contributes to the basioccipital bone and that the location of the *Pax7* lineage correlates with the morphology most effected by *Pax7* deficiency. Our results suggest that the *Pax7*-deficient mouse may be a suitable model for investigating the genetic control of the location and orientation of the foramen magnum, and changes in the breadth of the basioccipital.

M. Chen, G. Grinstein, C. R. Johnson, J. Kennedy, M. Tory.
**“Pathways for Theoretical Advances in Visualization,”** In *IEEE Computer Graphics and Applications*, IEEE, pp. 103--112. July, 2017.

More than a decade ago, Chris Johnson proposed the "Theory of Visualization" as one of the top research problems in visualization. Since then, there have been several theory-focused events, including three workshops and three panels at IEEE Visualization (VIS) Conferences. Together, these events have produced a set of convincing arguments.

M. Feiszli, A. Narayan.
**“Numerical Computation of Weil-Peterson Geodesics in the Universal Teichmueller Space,”** In *SIAM Journal on Imaging Sciences*, Vol. 10, No. 3, SIAM, pp. 1322--1345. Jan, 2017.

DOI: 10.1137/15M1043947

We propose an optimization algorithm for computing geodesics on the universal Teichm\"uller space T(1) in the Weil-Petersson (WP) metric. Another realization for T(1) is the space of planar shapes, modulo translation and scale, and thus our algorithm addresses a fundamental problem in computer vision: compute the distance between two given shapes. The identification of smooth shapes with elements on T(1) allows us to represent a shape as a diffeomorphism on S1. Then given two diffeomorphisms on S1 (i.e., two shapes we want connect with a flow), we formulate a discretized WP energy and the resulting problem is a boundary-value minimization problem. We numerically solve this problem, providing several examples of geodesic flow on the space of shapes, and verifying mathematical properties of T(1). Our algorithm is more general than the application here in the sense that it can be used to compute geodesics on any other Riemannian manifold.

C. Gritton, J. Guilkey, J. Hooper, D. Bedrov, R. M. Kirby, M. Berzins.
**“Using the material point method to model chemical/mechanical coupling in the deformation of a silicon anode,”** In *Modelling and Simulation in Materials Science and Engineering*, Vol. 25, No. 4, pp. 045005. 2017.

The lithiation and delithiation of a silicon battery anode is modeled using the material point method (MPM). The main challenges in modeling this process using the MPM is to simulate stress dependent diffusion coupled with concentration dependent stress within a material that undergoes large deformations. MPM is chosen as the numerical method of choice because of its ability to handle large deformations. A method for modeling diffusion within MPM is described. A stress dependent model for diffusivity and three different constitutive models that fully couple the equations for stress with the equations for diffusion are considered. Verifications tests for the accuracy of the numerical implementations of the models and validation tests with experimental results show the accuracy of the approach. The application of the fully coupled stress diffusion model implemented in MPM is applied to modeling the lithiation and delithiation of silicon nanopillars.

L. Guo, A. Narayan, T. Zhou, Y. Chen.
**“Stochastic Collocation Methods via L1 Minimization Using Randomized Quadratures,”** In *SIAM Journal on Scientific Computing*, Vol. 39, No. 1, pp. A333--A359. Jan, 2017.

ISSN: 1064-8275

DOI: 10.1137/16M1059680

In this work, we discuss the problem of approximating a multivariate function via ℓ1 minimization method, using a random chosen sub-grid of the corresponding tensor grid of Gaussian points. The independent variables of the function are assumed to be random variables, and thus, the framework provides a non-intrusive way to construct the generalized polynomial chaos expansions, stemming from the motivating application of Uncertainty Quantification (UQ). We provide theoretical analysis on the validity of the approach. The framework includes both the bounded measures such as the uniform and the Chebyshev measure, and the unbounded measures which include the Gaussian measure. Several numerical examples are given to confirm the theoretical results.

J. K. Holmen, A. Humphrey, D. Sutherland, M. Berzins.
**“Improving Uintah's Scalability Through the Use of Portable Kokkos-Based Data Parallel Tasks,”** In *Proceedings of the Practice and Experience in Advanced Research Computing 2017 on Sustainability, Success and Impact*, PEARC17, No. 27, pp. 27:1--27:8. 2017.

ISBN: 978-1-4503-5272-7

DOI: 10.1145/3093338.3093388

The University of Utah's Carbon Capture Multidisciplinary Simulation Center (CCMSC) is using the Uintah Computational Framework to predict performance of a 1000 MWe ultra-supercritical clean coal boiler. The center aims to utilize the Intel Xeon Phi-based DOE systems, Theta and Aurora, through the Aurora Early Science Program by using the Kokkos C++ library to enable node-level performance portability. This paper describes infrastructure advancements and portability improvements made possible by our integration of Kokkos within Uintah. Scalability results are presented that compare serial and data parallel task execution models for a challenging radiative heat transfer calculation, central to the center's predictive boiler simulations. These results demonstrate both good strong-scaling characteristics to 256 Knights Landing (KNL) processors on the NSF Stampede system, and show the KNL-based calculation to compete with prior GPU-based results for the same calculation.

J. Jakeman, A. Narayan, T. Zhou.
**“A Generalized Sampling and Preconditioning Scheme for Sparse Approximation of Polynomial Chaos Expansions,”** In *SIAM Journal on Scientific Computing*, Vol. 39, No. 3, SIAM, pp. A1114--A1144. Jan, 2017.

ISSN: 1064-8275

DOI: 10.1137/16M1063885

In this paper we propose an algorithm for recovering sparse orthogonal polynomials using stochastic collocation. Our approach is motivated by the desire to use generalized polynomial chaos expansions (PCE) to quantify uncertainty in models subject to uncertain input parameters. The standard sampling approach for recovering sparse polynomials is to use Monte Carlo (MC) sampling of the density of orthogonality. However MC methods result in poor function recovery when the polynomial degree is high. Here we propose a general algorithm that can be applied to any admissible weight function on a bounded domain and a wide class of exponential weight functions defined on unbounded domains. Our proposed algorithm samples with respect to the weighted equilibrium measure of the parametric domain, and subsequently solves a preconditioned ℓ1-minimization problem, where the weights of the diagonal preconditioning matrix are given by evaluations of the Christoffel function. We present theoretical analysis to motivate the algorithm, and numerical results that show our method is superior to standard Monte Carlo methods in many situations of interest. Numerical examples are also provided that demonstrate that our proposed Christoffel Sparse Approximation algorithm leads to comparable or improved accuracy even when compared with Legendre and Hermite specific algorithms.

J. Jiang, Y. Chen, A. Narayan.
**“Offline-Enhanced Reduced Basis Method through adaptive construction of the Surrogate Parameter Domain,”** In *Journal of Scientific Computing*, Vol. 73, No. 2-3, pp. 853--875. 2017.

ISSN: 0885-7474, 1573-7691

DOI: 10.1007/s10915-017-0551-3

The Reduced Basis Method (RBM) is a popular certified model reduction approach for solving parametrized partial differential equations. One critical stage of the \textitoffline portion of the algorithm is a greedy algorithm, requiring maximization of an error estimate over parameter space. In practice this maximization is usually performed by replacing the parameter domain continuum with a discrete "training" set. When the dimension of parameter space is large, it is necessary to significantly increase the size of this training set in order to effectively search parameter space. Large training sets diminish the attractiveness of RBM algorithms since this proportionally increases the cost of the offline phase.

In this work we propose novel strategies for offline RBM algorithms that mitigate the computational difficulty of maximizing error estimates over a training set. The main idea is to identify a subset of the training set, a "surrogate parameter domain" (SPD), on which to perform greedy algorithms. The SPD's we construct are much smaller in size than the full training set, yet our examples suggest that they are accurate enough to represent the solution manifold of interest at the current offline RBM iteration. We propose two algorithms to construct the SPD: Our first algorithm, the Successive Maximization Method (SMM) method, is inspired by inverse transform sampling for non-standard univariate probability distributions. The second constructs an SPD by identifying pivots in the Cholesky Decomposition of an approximate error correlation matrix. We demonstrate the algorithm through numerical experiments, showing that the algorithm is capable of accelerating offline RBM procedures without degrading accuracy, assuming that the solution manifold has low Kolmogorov width.

M. Kern, A. Lex, N. Gehlenborg, C. R. Johnson.
**“Interactive Visual Exploration And Refinement Of Cluster Assignments,”** In *BMC Bioinformatics*, Cold Spring Harbor Laboratory, April, 2017.

DOI: 10.1101/123844

**Background:**

With ever-increasing amounts of data produced in biology research, scientists are in need of efficient data analysis methods. Cluster analysis, combined with visualization of the results, is one such method that can be used to make sense of large data volumes. At the same time, cluster analysis is known to be imperfect and depends on the choice of algorithms, parameters, and distance measures. Most clustering algorithms don't properly account for ambiguity in the source data, as records are often assigned to discrete clusters, even if an assignment is unclear. While there are metrics and visualization techniques that allow analysts to compare clusterings or to judge cluster quality, there is no comprehensive method that allows analysts to evaluate, compare, and refine cluster assignments based on the source data, derived scores, and contextual data.**Results:**

In this paper, we introduce a method that explicitly visualizes the quality of cluster assignments, allows comparisons of clustering results and enables analysts to manually curate and refine cluster assignments. Our methods are applicable to matrix data clustered with partitional, hierarchical, and fuzzy clustering algorithms. Furthermore, we enable analysts to explore clustering results in context of other data, for example, to observe whether a clustering of genomic data results in a meaningful differentiation in phenotypes.**Conclusions:**

Our methods are integrated into Caleydo StratomeX, a popular, web-based, disease subtype analysis tool. We show in a usage scenario that our approach can reveal ambiguities in cluster assignments and produce improved clusterings that better differentiate genotypes and phenotypes.

A. Narayan, J. Jakeman, T. Zhou.
**“A Christoffel function weighted least squares algorithm for collocation approximations,”** In *Mathematics of Computation*, Vol. 86, No. 306, pp. 1913--1947. 2017.

ISSN: 0025-5718, 1088-6842

DOI: 10.1090/mcom/3192

We propose, theoretically investigate, and numerically validate an algorithm for the Monte Carlo solution of least-squares polynomial approximation problems in a collocation frame- work. Our method is motivated by generalized Polynomial Chaos approximation in uncertainty quantification where a polynomial approximation is formed from a combination of orthogonal polynomials. A standard Monte Carlo approach would draw samples according to the density of orthogonality. Our proposed algorithm samples with respect to the equilibrium measure of the parametric domain, and subsequently solves a weighted least-squares problem, with weights given by evaluations of the Christoffel function. We present theoretical analysis to motivate the algorithm, and numerical results that show our method is superior to standard Monte Carlo methods in many situations of interest.

T.A.J. Ouermi, A. Knoll, R.M. Kirby, M. Berzins.
**“OpenMP 4 Fortran Modernization of WSM6 for KNL,”** In *Proceedings of the Practice and Experience in Advanced Research Computing 2017 on Sustainability, Success and Impact*, PEARC17, No. 12, ACM, pp. 12:1--12:8. 2017.

ISBN: 978-1-4503-5272-7

DOI: 10.1145/3093338.3093387

Parallel code portability in the petascale era requires modifying existing codes to support new architectures with large core counts and SIMD vector units. OpenMP is a well established and increasingly supported vehicle for portable parallelization. As architectures mature and compiler OpenMP implementations evolve, best practices for code modernization change as well. In this paper, we examine the impact of newer OpenMP features (in particular OMP SIMD) on the Intel Xeon Phi Knights Landing (KNL) architecture, applied in optimizing loops in the single moment 6-class microphysics module (WSM6) in the US Navy's NEPTUNE code. We find that with functioning OMP SIMD constructs, low thread invocation overhead on KNL and reduced penalty for unaligned access compared to previous architectures, one can leverage OpenMP 4 to achieve reasonable scalability with relatively minor reorganization of a production physics code.

T.A.J. Ouermi, A. Knoll, R.M. Kirby, M. Berzins.
**“Optimization Strategies for WRF Single-Moment 6-Class Microphysics Scheme (WSM6) on Intel Microarchitectures,”** In *Proceedings of the fifth international symposium on computing and networking (CANDAR 17). Awarded Best Paper *, IEEE, 2017.

Optimizations in the petascale era require modifications of existing codes to take advantage of new architectures with large core counts and SIMD vector units. This paper examines high-level and low-level optimization strategies for numerical weather prediction (NWP) codes. These strategies employ thread-local structures of arrays (SOA) and an OpenMP directive such as OMP SIMD. These optimization approaches are applied to the Weather Research Forecasting single-moment 6-class microphysics schemes (WSM6) in the US Navy NEPTUNE system. The results of this study indicate that the high-level approach with SOA and low-level OMP SIMD improves thread and vector parallelism by increasing data and temporal locality. The modified version of WSM6 runs 70x faster than the original serial code. This improvement is about 23.3x faster than the performance achieved by Ouermi et al., and 14.9x faster than the performance achieved by Michalakes et al.

B. Peterson, A. Humphrey, J. Schmidt, M. Berzins.
**“Addressing Global Data Dependencies in Heterogeneous Asynchronous Runtime Systems on GPUs. Awarded Best Paper,”** In *Proceedings of the Third International Workshop on Extreme Scale Programming Models and Middleware - ESPM2'17*, ACM, 2017.

DOI: 10.1145/3152041.3152082

Large-scale parallel applications with complex global data dependencies beyond those of reductions pose significant scalability challenges in an asynchronous runtime system. Internodal challenges include identifying the all-to-all communication of data dependencies among the nodes. Intranodal challenges include gathering together these data dependencies into usable data objects while avoiding data duplication. This paper addresses these challenges within the context of a large-scale, industrial coal boiler simulation using the Uintah asynchronous many-task runtime system on GPU architectures. We show significant reduction in time spent analyzing data dependencies through refinements in our dependency search algorithm. Multiple task graphs are used to eliminate subsequent analysis when task graphs change in predictable and repeatable ways. Using a combined data store and task scheduler redesign reduces data dependency duplication ensuring that problems fit within host and GPU memory. These modifications did not require any changes to application code or sweeping changes to the Uintah runtime system. We report results running on the DOE Titan system on 119K CPU cores and 7.5K GPUs simultaneously. Our solutions can be generalized to other task dependency problems with global dependencies among thousands of nodes which must be processed efficiently at large scale.

P. Seshadri, A. Narayan, S. Mahadevan.
**“Effectively Subsampled Quadratures for Least Squares Polynomial Approximations,”** In *SIAM/ASA Journal on Uncertainty Quantification*, pp. 1003--1023. Jan, 2017.

This paper proposes a new deterministic sampling strategy for constructing polynomial chaos approximations for expensive physics simulation models. The proposed approach, effectively subsampled quadratures involves sparsely subsampling an existing tensor grid using QR column pivoting. For polynomial interpolation using hyperbolic or total order sets, we then solve the following square least squares problem. For polynomial approximation, we use a column pruning heuristic that removes columns based on the highest total orders and then solves the tall least squares problem. While we provide bounds on the condition number of such tall submatrices, it is difficult to ascertain how column pruning effects solution accuracy as this is problem specific. We conclude with numerical experiments on an analytical function and a model piston problem that show the efficacy of our approach compared with randomized subsampling. We also show an example where this method fails.