Search

We consider numerical simulation of the isotropic elastic wave equations arising from seismic applications with non-trivial land topography. The more flexible finite element method is applied to the shallow region of the simulation domain to account for the topography, and combined with the more efficient finite difference method that is applied to the deep region of the simulation domain. We demonstrate that these two discretization methods, albeit starting from different formulations of the elastic wave equation, can be joined together smoothly via weakly imposed interface conditions. Discrete energy analysis is employed to derive the proper interface treatment, leading to an overall discretization that is energy-conserving. Numerical examples are presented to demonstrate the efficacy of the proposed interface treatment.

We present a robust and scalable preconditioner for the solution of large-scale linear systems that arise from the discretization of elliptic PDEs amenable to rank compression. The preconditioner is based on hierarchical low-rank approximations and the cyclic reduction method. The setup and application phases of the preconditioner achieve log-linear complexity in memory footprint and number of operations, and numerical experiments exhibit good weak and strong scalability at large processor counts in a distributed memory environment. Numerical experiments with linear systems that feature symmetry and nonsymmetry, definiteness and indefiniteness, constant and variable coefficients demonstrate the preconditioner applicability and robustness. Furthermore, it is possible to control the number of iterations via the accuracy threshold of the hierarchical matrix approximations and their arithmetic operations, and the tuning of the admissibility condition parameter. Together, these parameters allow for optimization of the memory requirements and performance of the preconditioner., Comment: 24 pages, Elsevier Journal of Computational and Applied Mathematics, Dec 2017

The sum of Kronecker products (SKP) representation for spatial covariance matrices from gridded observations and a corresponding adaptive-cross-approximation-based framework for building the Kronecker factors are investigated. The time cost for constructing an n -dimensional covariance matrix is O ( n k 2 ) and the total memory footprint is O ( n k ) , where k is the number of Kronecker factors. The memory footprint under the SKP representation is compared with that under the hierarchical representation and found to be one order of magnitude smaller. A Cholesky factorization algorithm under the SKP representation is proposed and shown to factorize a one-million dimensional covariance matrix in under 600 seconds on a standard scientific workstation. With the computed Cholesky factor, simulations of Gaussian random fields in one million dimensions can be achieved at a low cost for a wide range of spatial covariance functions.

Shared memory parallelization of the flux kernel of PETSc-FUN3D, an unstructured tetrahedral mesh Euler flow code previously studied for distributed memory and multi-core shared memory, is evaluated on up to 61 cores per node and up to 4 threads per core. We explore several thread-level optimizations to improve flux kernel performance on the state-of-the-art many integrated core (MIC) Intel processor Xeon Phi “Knights Corner,” with a focus on strong thread scaling. While the linear algebraic kernel is bottlenecked by memory bandwidth for even modest numbers of cores sharing a common memory, the flux kernel, which arises in the control volume discretization of the conservation law residuals and in the formation of the preconditioner for the Jacobian by finite-differencing the conservation law residuals, is compute-intensive and is known to exploit effectively contemporary multi-core hardware. We extend study of the performance of the flux kernel to the Xeon Phi in three thread affinity modes, namely scatter, compact, and balanced, in both offload and native mode, with and without various code optimizations to improve alignment and reduce cache coherency penalties. Relative to baseline “out-of-the-box” optimized compilation, code restructuring optimizations provide about 3.8x speedup using the offload mode and about 5x speedup using the native mode. Even with these gains for the flux kernel, with respect to execution time the MIC simply achieves par with optimized compilation on a contemporary multi-core Intel CPU, the 16-core Sandy Bridge E5 2670. Nevertheless, the optimizations employed to reduce the data motion and cache coherency protocol penalties of the MIC are expected to be of value for CFD and many other unstructured applications as many-core architecture evolves. We explore large-scale distributed-shared memory performance on the Cray XC40 supercomputer, to demonstrate that optimizations employed on Phi hybridize to this context, where each of thousands of nodes are comprised of two sockets of Intel Xeon Haswell CPUs with 32 cores per node.

In computational fluid dynamics, the demand for increasingly multidisciplinary reliable simulations, for both analysis and design optimization purposes, requires transformational advances in individual components of future solvers. At the algorithmic level, hardware compatibility and efficiency are of paramount importance in determining viability at exascale and beyond. However, equally important (if not more so) is algorithmic robustness with minimal user intervention, which becomes progressively more challenging to achieve as problem size and physics complexity increase. We numerically show that low and high order entropy stable collocated discontinuous Galerkin discretizations based on summation-by-part operators and simultaneous-approximation-terms technique provide an essential step toward a truly enabling technology in terms of reliability and robustness for both under-resolved turbulent flow simulations and flows with discontinuities.

Simulation of subsurface flows in porous media is difficult due to the nonlinearity of the operators and the high heterogeneity of material coefficients. In this paper, we present a scalable fully implicit solver for incompressible two-phase flows based on overlapping domain decomposition methods. Specifically, an inexact Newton-Krylov algorithm with analytic Jacobian is used to solve the nonlinear systems arising from the discontinuous Galerkin discretization of the governing equations on 3D unstructured grids. The linear Jacobian system is preconditioned by additive Schwarz algorithms, which are naturally suitable for parallel computing. We propose a hybrid two-level version of the additive Schwarz preconditioner consisting of a nested coarse space to improve the robustness and scalability of the classical one-level version. On the coarse level, a smaller linear system arising from the same discretization of the problem on a coarse grid is solved by using GMRES with a one-level preconditioner until a relative tolerance is reached. Numerical experiments are presented to demonstrate the effectiveness and efficiency of the proposed solver for 3D heterogeneous medium problems. We also report the parallel scalability of the proposed algorithms on a supercomputer with up to 8 , 192 processor cores.

To further enhance the capacity of next generation wireless communication systems, massive multiple-input multiple-output (MIMO) has recently appeared as a necessary enabling technology to achieve high performance signal processing for large-scale multiple antennas. However, massive MIMO systems inevitably generate signal processing overheads, which translate into ever-increasing rate of complexity, and therefore, such systems may not maintain the inherent real-time requirement of wireless systems. We redesign the nonlinear sphere decoder method to increase the performance of the system, cast most memory-bound computations into compute-bound operations to reduce the overall complexity, and maintain the real-time processing thanks to the graphics processing unit (GPU) computational power. We show a comprehensive complexity and performance analysis on an unprecedented MIMO system scale, which can ease the design phase toward simulating future massive MIMO wireless systems.

We present a numerical algorithm for simulating the spinodal decomposition described by the three dimensional Cahn-Hilliard-Cook (CHC) equation, which is a fourth-order stochastic partial differential equation with a noise term. The equation is discretized in space and time based on a fully implicit, cell-centered finite difference scheme, with an adaptive time-stepping strategy designed to accelerate the progress to equilibrium. At each time step, a parallel Newton-Krylov-Schwarz algorithm is used to solve the nonlinear system. We discuss various numerical and computational challenges associated with the method. The numerical scheme is validated by a comparison with an explicit scheme of high accuracy (and unreasonably high cost). We present steady state solutions of the CHC equation in two and three dimensions. The effect of the thermal fluctuation on the spinodal decomposition process is studied. We show that the existence of the thermal fluctuation accelerates the spinodal decomposition process and that the final steady morphology is sensitive to the stochastic noise. We also show the evolution of the energies and statistical moments. In terms of the parallel performance, it is found that the implicit domain decomposition approach scales well on supercomputers with a large number of processors.

We study a new optimization scheme that generates smooth and robust solutions for Dirichlet velocity boundary control (DVBC) of conjugate heat transfer (CHT) processes. The solutions to the DVBC of the incompressible Navier-Stokes equations are typically nonsmooth, due to the regularity degradation of the boundary stress in the adjoint Navier-Stokes equations. This nonsmoothness is inherited by the solutions to the DVBC of CHT processes, since the CHT process couples the Navier-Stokes equations of fluid motion with the convection-diffusion equations of fluid-solid thermal interaction. Our objective in the CHT boundary control problem is to select optimally the fluid inflow profile that minimizes an objective function that involves the sum of the mismatch between the temperature distribution in the fluid system and a prescribed temperature profile and the cost of the control.Our strategy to resolve the nonsmoothness of the boundary control solution is based on two features, namely, the objective function with a regularization term on the gradient of the control profile on both the continuous and the discrete levels, and the optimization scheme with either explicit or implicit smoothing effects, such as the smoothed Steepest Descent and the Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) methods. Our strategy to achieve the robustness of the solution process is based on combining the smoothed optimization scheme with the numerical continuation technique on the regularization parameters in the objective function. In the section of numerical studies, we present two suites of experiments. In the first one, we demonstrate the feasibility and effectiveness of our numerical schemes in recovering the boundary control profile of the standard case of a Poiseuille flow. In the second one, we illustrate the robustness of our optimization schemes via solving more challenging DVBC problems for both the channel flow and the flow past a square cylinder, which use initial control profiles far from optimal and require the numerical continuation technique applied on regularization parameters. We believe our solution strategy is general and can be applied to other large-scale optimal control problems which involve multiphysics processes and require smooth approximations to the optimal control profile.

Numerical simulations of the four-field extended magnetohydrodynamics (MHD) equations with hyper-resistivity terms present a difficult challenge because of demanding spatial resolution requirements. A time-dependent sequence of r-refinement adaptive grids obtained from solving a single Monge-Ampere (MA) equation addresses the high-resolution requirements near the x-point for numerical simulation of the magnetic reconnection problem. The MHD equations are transformed from Cartesian coordinates to solution-defined curvilinear coordinates. After the application of an implicit scheme to the time-dependent problem, the parallel Newton-Krylov-Schwarz (NKS) algorithm is used to solve the system at each time step. Convergence and accuracy studies show that the curvilinear solution requires less computational effort than a pure Cartesian treatment. This is due both to the more optimal placement of the grid points and to the improved convergence of the implicit solver, nonlinearly and linearly. The latter effect, which is significant (more than an order of magnitude in number of inner linear iterations for equivalent accuracy), does not yet seem to be widely appreciated.

The best paths to the exascale summit are debatable, but all are narrow and treacherous, constrained by fundamental laws of physics, capital cost, operating cost, power requirements, programmability, and reliability. Many scientific and engineering applications force the modeling community to attempt to scale this summit. Drawing on vendor projections and experiences with scientific codes on contemporary platforms, we outline the challenges and propose roles and essential adaptations for mathematical modelers in one of the great global scientific quests the next decade.

This paper presents a set of computationally efficient, adaptive grids for magnetic reconnection phenomenon where the current density can develop large gradients in the reconnection region. Four-field extended MagnetoHydroDynamics (MHD) equations with hyperviscosity terms are transformed so that the curvilinear coordinates replace the Cartesian coordinates as the independent variables, and moving grids' velocities are also considered in this transformed system as a part of interpolating the physical solutions from the old grid to the new grid as time advances. The curvilinear coordinates derived from the current density through the Monge–Kantorovich (MK) optimization approach help to reduce the resolution requirements during the computation.

A parallel, fully coupled, nonlinearly implicit Newton-Krylov-Schwarz algorithm is proposed for the numerical simulation of a magnetic reconnection problem described by a system of resistive Hall magnetohydrodynamics equations in slab geometry. A key component of the algorithm is a restricted additive Schwarz preconditioner defined for problems with doubly periodic boundary conditions. We show numerically that with such a preconditioned nonlinearly implicit method the time step size is no longer constrained by the CFL number or the convergence of the Newton solver. We report the parallel performance of the algorithm and software on machines with thousands of processors.

Jacobian-free Newton-Krylov (JFNK) methods are synergistic combinations of Newton-type methods for superlinearly convergent solution of nonlinear equations and Krylov subspace methods for solving the Newton correction equations. The link between the two methods is the Jacobian-vector product, which may be probed approximately without forming and storing the elements of the true Jacobian, through a variety of means. Various approximations to the Jacobian matrix may still be required for preconditioning the resulting Krylov iteration. As with Krylov methods for linear problems, successful application of the JFNK method to any given problem is dependent on adequate preconditioning. JFNK has potential for application throughout problems governed by nonlinear partial differential equations and integro-differential equations. In this survey paper, we place JFNK in context with other nonlinear solution algorithms for both boundary value problems (BVPs) and initial value problems (IVPs). We provide an overview of the mechanics of JFNK and attempt to illustrate the wide variety of preconditioning options available. It is emphasized that JFNK can be wrapped (as an accelerator) around another nonlinear fixed point method (interpreted as a preconditioning process, potentially with significant code reuse). The aim of this paper is not to trace fully the evolution of JFNK, nor to provide proofs of accuracy or optimal convergence for all of the constituent methods, but rather to present the reader with a perspective on how JFNK may be applicable to applications of interest and to provide sources of further practical information.

Fluid dynamical simulations based on finite discretizations on (quasi-)static grids scale well in parallel, but execute at a disappointing percentage of per-processor peak floating point operation rates without special attention to layout and access ordering of data. We document both claims from our experience with an unstructured grid CFD code that is typical of the state of the practice at NASA. These basic performance characteristics of PDE-based codes can be understood with surprisingly simple models, for which we quote earlier work, presenting primarily experimental results. The performance models and experimental results motivate algorithmic and software practices that lead to improvements in both parallel scalability and per node performance. This snapshot of ongoing work updates our 1999 Bell Prize-winning simulation on ASCI computers.

Two paradigms for distributed-memory parallel computation that free the application programmer from the details of message passing are compared for an archetypal structured scientific computation --- a nonlinear, structured-grid partial differential equation boundary value problem --- using the same algorithm on the same hardware. Both paradigms, parallel libraries represented by Argonne''s PETSc, and parallel languages represented by the Portland Group''s HPF, are found to be easy to use for this problem class, and both are reasonably effective in exploiting concurrency after a short learning curve. The level of involvement required by the application programmer under either paradigm includes specification of the data partitioning (corresponding to a geometrically simple decomposition of the domain of the PDE). Programming in SPMD style for the PETSc library requires writing the routines that discretize the PDE and its Jacobian, managing subdomain-to-processor mappings (affine global-to-local index mappings), and interfacing to library solver routines. Programming for HPF requires a complete sequential implementation of the same algorithm, introducing concurrency through subdomain blocking (an effort similar to the index mapping), and modest experimentation with rewriting loops to elucidate to the compiler the latent concurrency. Correctness and scalability are cross-validated on up to 32 nodes of an IBM SP2.

We introduce a model for communication costs in parallel processing environments, called the hyperbolic model, which generalizes two-parameter dedicated-link models in an analytically simple way. The communication system is modeled as a directed communication graph in which terminal nodes represent the application processes and internal nodes, called communication blocks (CBs), reflect the layered structure of the underlying communication architecture. ACBis characterized by a two-parameter hyperbolic function of the message size that represents the service time needed for processing the message. Rules are given for reducing a communication graph consisting of manyCBs to an equivalent two-parameter form, while maintaining a good approximation for the service time. We demonstrate a tight fit of the estimates of the cost of communication based on the model with actual measurements of the communication and synchronization time between end processes. We compare the hyperbolic model with other two-parameter models and, in appropriate limits, show its compatibility with the LogP model.

Nonlinear elastoplastic boundary element analysis (BEA) involves an algebraic subproblem requiring the solution of dense nonsymmetric matrix equations with an evolving right hand side vector. When multiple right hand side vectors are present, direct matrix triangular factorization techniques have been the compelling choice, amortizing the work of a single matrix factorization over the sequence of multiple fast ‘solutions’ of the resulting triangular systems. Recently, the superior performance of preconditioned Krylov equation solvers in linear BEA has also been documented. In this paper, the superior performance of preconditioned Krylov equation solvers is shown to be extendable to elastoplastic BEA. This is accomplished by exploiting the strategic reuse of the preconditioner, its factorization, and the Krylov vectors computed in the solution for the first right hand side vector, in the subsequent solution of matrix equations with multiple ‘nearby’ right hand side vectors. The details associated with this strategy are given, and the computer resources required in three dimensional elastoplastic BEA are used to quantify the computational efficiency associated with this new algorithm.

Parallel block-preconditioned domain-decomposed Krylov methods for sparse linear systems are described and illustrated on large two-dimensional model problems and Jacobian matrices from different stages of a nonlinear multicomponent problem in chemically reacting flows. The main motivation of the work is to examine the practicality of parallelization, under the domain decomposition paradigm, of the solution of systems of equations typical of implicit finite difference applications from fluid dynamics. Such systems presently lie beyond the realm of most of the theory for domain-decomposed symmetric or nonsymmetric scalar operators. We describe techniques depending formally only on the sparsity structure of the linear operator and thus of broad applicability. Results of tests run on an Encore Multimax with up to 16 processors demonstrate their utility in the coarse-granularity parallelization of hydrocodes: parallel efficiencies in the 40 to 100 percent range are available on the largest number of processors employed over a mix of problems, relative to a serial approach employing the same iterative technique (GMRES) and preconditioner (ILU) on a single domain. These efficiencies are already competitive with results from undecomposed parallel implementations of ILU-preconditioned GMRES on the same multiprocessor, and many avenues for their improvement remain unexplored.

We consider the application of domain decomposition techniques to the solution of sparse linear systems arising from implicit PDE discretizations on parallel computers. Representatives of two popular MIMD architectures, message passing (the Intel iPSC/2-SX) and shared memory (the Encore Multimax 320), are employed. We run the same numerical experiments on each, namely stripwise and boxwise decompositions of the unit square, using up to 64 subdomains and containing up to 64K degrees of freedom. We produce a tight-fitting complexity model for the former and discuss the difficulty of doing so for the latter. We also evaluate which of three types of domain decomposition preconditioners that have appeared in the literature of self-adjoint elliptic problems are most efficient in different regions of machine-problem parameter space. Some form of global sharing of information in the preconditioner is required for efficient overall parallel implementation in the region of most practical interest (large problem sizes and large numbers of processors); otherwise, an increasing iteration count inveighs against the gains of concurrency. Our results on a per iteration basis also hold for sparse discrete systems arising from other types of partial differential equations, but in the absence of a theory for the dependence of the convergence rate upon the granularity of the decomposition, the overall results are only suggestive for more general systems.