SUNDIALS Multiphysics+MPIManyVector Performance Testing
Daniel R. Reynolds David J. Gardner Cody J. Balos Carol S. Woodward


empty




empty








Generating Fair Universal Representations using Adversarial Models
    Yusuf Said Eroglu, Fatih Erden, and Ismail Guvenc

Department of Electrical and Computer Engineering, North Carolina State University, Raleigh, NC

Email: {yeroglu, ferden, iguvenc}@ncsu.edu
This work is supported in part by NSF CNS award 1422062.

    December 30, 2023
============================================================================================================================================================================================================================================================





§ INTRODUCTION



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 <cit.> 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:


  * SUNDIALS' new  module, that allows extreme
  flexibility in how a solution “vector” is staged on computational
  resources.

  * ARKode's new multirate integration module, ,
  allowing high-order accurate calculations that subcycle “fast”
  processes within “slow” ones. 

  * SUNDIALS' new flexible linear solver interfaces, that allow
  streamlined specification of problem-specific linear solvers. 

  * SUNDIALS' new  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 
   call). 

We anticipate that subsequent reports will extend this work to
investigate a variety of other new features, including SUNDIALS' 
generic  interface and accelerator-enabled 
 modules, and upcoming 
extensions to support custom “fast” integrators (that leverage
problem structure) and IMEX integration of the “slow” time scale (to
add diffusion).




§ PROBLEM DESCRIPTION



We simulate the three-dimensional nonlinear inviscid compressible Euler
equations, combined with advection and reaction of chemical species,

    _t = -∇·() + () + (,t).

Here, the independent variables are (,t) = (x,y,z,t) ∈Ω× (t_0,
t_f], where the spatial domain is a three-dimensional cube, Ω =
[x_l, x_r] × [y_l, y_r] × [z_l,z_r].  The partial differential
equation is completed using initial condition (,t_0) = _0() and
face-specific boundary conditions, {xlbc,  xrbc,  ylbc,  yrbc,  zlbc, 
zrbc}, corresponding to conditions that may be separately applied at
the spatial locations {(x_l,y,z), (x_r,y,z), (x,y_l,z),
(x,y_r,z), (x,y,z_l), (x,y,z_r)}, respectively, and where each
condition may be any one of:


  * periodic (requires that both faces in this direction use
  this condition),

  * homogeneous Neumann (i.e., ∇ w_i(,t)·n=0
  for ∈∂Ω with outward-normal vector 𝐧, 
  and for each species w_i), 

  * homogeneous Dirichlet (i.e., w_i(,t) = 0 for
  ∈∂Ω and for each species w_i), or 

  * reflecting (i.e., homogeneous Neumann for all species
  except the momentum field perpendicular to that face, that
  has a homogeneous Dirichlet condition), 

The computed solution is given by
= [     ρ ρ v_x ρ v_y ρ v_z   e_t       ]^T
= [   ρ m_x m_y m_z e_t     ]^T,
that corresponds to the density (ρ), x,y,z-momentum
(m_x,m_y,m_z), total energy per unit volume (e_t), and vector of
chemical densities (∈^n_c) that are advected
along with the fluid.  The advective fluxes = [ F_x F_y F_z ]^T are given by 

    F_x()    = [       ρ v_x ρ v_x^2 + p   ρ v_x v_y   ρ v_x v_z v_x (e_t+p)         v_x ]^T
    
      F_y()    = [       ρ v_y   ρ v_x v_y ρ v_y^2 + p   ρ v_y v_z v_y (e_t+p)         v_y ]^T
    
      F_z()    = [       ρ v_z   ρ v_x v_z   ρ v_y v_z ρ v_z^2 + p v_z (e_t+p)         v_z ]^T.

The reaction term () and external force (,t) are
test-problem-dependent, and the ideal gas equation of state relates the pressure
and total energy density,

    p    = R/c_v(e_t - ρ/2 (v_x^2 + v_y^2 + v_z^2))
    ⇔       
    
      e_t    = p c_v/R + ρ/2(v_x^2 + v_y^2 + v_z^2),

or equivalently,

    p    = (γ-1) (e_t - ρ/2 (v_x^2 + v_y^2 + v_z^2))
    ⇔       
    
      e_t    = p/γ-1 + ρ/2(v_x^2 + v_y^2 + v_z^2),

The above model includes the physical parameters:


  * R is the specific ideal gas constant (287.14 J/kg/K for air).

  * c_v is the specific heat capacity at constant volume (717.5
  J/kg/K for air),

  * γ is the ratio of specific heats, γ = c_p/c_v =
  1 + R/c_v; this is typically 1.4 for air, and 5/3 for astrophysical
  gases.

The speed of sound in the gas is given by

    c = √(γ p/ρ).

The fluid variables (ρ, 𝐦, and e_t) are non-dimensionalized; 
when converted to physical CGS values these have units:


  * [ρ] = g/cm^3,

  * [v_x] = [v_y] = [v_z] = cm/s, which implies that 
  [m_x] = [m_y] = [m_z] = g/cm^2/s,

  * [e_t] = g/cm/s^2.

The chemical densities have physical units [_i] = g/cm^3, although when
these are transported by the fluid these are converted to dimensionless units as well.





§ IMPLEMENTATION



We apply a method of lines approach for converting the system of partial 
differential equations (PDEs) (<ref>) into a discrete set 
of equations.  To this end, we first discretize in space, converting the 
PDE system into a very large system of ordinary differential equation 
(ODE) initial-value problems (IVPs).  We then apply the  
time-integrator from the ARKode SUNDIALS package.  This time 
integration approach in turn requires a sub-integrator for the “fast” 
chemical reactions, for which we employ the  
time-integrator, also from ARKode.  , in turn, 
requires the solution of very large-scale systems of nonlinear 
algebraic equations.  We discuss our use of each of the above components
in the following subsections.




 §.§ Spatial discretization



We discretize the domain Ω into a uniform grid of
dimensions n_x × n_y × n_z, such that we have a three-dimensional
rectangular cuboid of cell-centered values (x_i,y_j,z_k) wherein

    x_i = x_l + (i + 1/2) Δ x,        Δ x = x_r-x_l/n_x,      i=0,…,n_x-1
    
      y_j = y_l + (j + 1/2) Δ y,        Δ y = y_r-y_l/n_y,      j=0,…,n_y-1,
    
      z_k = z_l + (k + 1/2) Δ z,        Δ z = z_r-z_l/n_z,      k=0,…,n_z-1.

This spatial domain is then decomposed in parallel using a standard 3D
domain decomposition approach over n_p MPI tasks, with layout
n_px× n_py× n_pz, defined automatically via the
 utility routine, as illustrated in Figure
<ref>.  This results in each MPI task “owning” a
local grid of dimensions n_xloc× n_yloc× n_zloc.



Within each cell in the domain we store 5+n_c variables,
corresponding to the values of  at that spatial location.  Here,
we employ the newly-introduced 
implementation, that allows creation of a single 
out of any valid set of subsidiary  objects.  To this
end, we store each of the five fluid fields (ρ, m_x, m_y,
m_z and e_t) in its own  object to
simplify access and I/O.  We then store all chemical species owned by each
MPI task in a single  object
(note: this will eventually be changed to use a device-specific
   object, such as ,
  , or ).  Each
MPI task then combines its pointers for the five fluid vectors, along
with its own chemical species vector, into its full “solution”
, , using the
 routine.  An illustration of this
 structure is shown in Figure <ref>.





The final item to note with regard to the spatial discretization is
our approach for approximating the flux divergence ∇·()
shown in equation (<ref>).  For this, we apply a
5th-order FD-WENO reconstruction, where we precisely follow the
algorithm laid out in the seminal paper by Shu <cit.>, which
we briefly summarize here.  In order to properly conserve mass,
momentum, and energy in (<ref>), we first compute the
fluxes at each of the six faces surrounding the cell (x_i,y_j,z_k),
and apply these in a standard conservative fashion, namely

    ∇·((x_i,y_j,z_k,t))
      ≈   1/Δ x[ F_x((x_i+1/2,y_j,z_k,t)) - F_x((x_i-1/2,y_j,z_k,t))] + 
       1/Δ y[ F_y((x_i,y_j+1/2,z_k,t)) - F_x((x_i,y_j-1/2,z_k,t))] + 
       1/Δ z[ F_z((x_i,y_j,z_k+1/2,t)) - F_x((x_i,y_j,z_k-1/2,t))].

Since the cells (x_i,y_j,z_k) and (x_i+1,y_j,z_k) share the same
flux value F_x((x_i+1/2,y_j,z_k,t)), this results in
conservation to full machine precision (modulo boundary conditions,
source terms, and reaction processes).  The 5th-order FD-WENO scheme is
used to construct each of these face-centered flux values.  This algorithm 
computes the flux F_x((x_i+1/2,y_j,z_k,t)) using a 6-point stencil of 
solution values (x_i-2,y_j,z_k,t), (x_i-1,y_j,z_k,t),
(x_i,y_j,z_k,t), (x_i+1,y_j,z_k,t),
(x_i+2,y_j,z_k,t), and (x_i+3,y_j,z_k,t), i.e., the
6 “closest” cell centers to the face (x_i+1/2,y_j,z_k) along the
x-direction.  The stencils are analagous in the y and z directions.  
Thus under our three-dimensional domain decomposition approach
outlined above, each MPI task must obtain three layers of “ghost” cells from
each of its six neighboring subdomains.

We thus compute the flux divergence ∇·() using the
following steps:


  * Begin exchange of boundary layers with neighbors via
  asynchronous  and  calls.

  * Compute and store the fluxes at each face in the strict interior
  of the subdomain, i.e., for all cell faces whose 6-point stencil
  involves no data from neighboring MPI tasks.  For example, we compute
  F_x((x_i-1/2,y_j,z_k)) over the ranges i=3,…,n_xloc-3.
  Since the vast majority of computational effort in this portion of the
  algorithm lies within the arithmetically intense FD-WENO 
  reconstruction itself, we first copy the local stencil of 
  6 × (5 + n_c) unknowns from their various  
  locations into a contiguous buffer before performing the FD-WENO
  reconstruction at each face.
  

  * Wait for completion of all asynchronous  and
   calls.

  * Compute the face-centered fluxes near subdomain boundaries, i.e., 
  those face-valued fluxes that were omitted above.  As these stencils now 
  depend on values from neighboring subdomains, the only difference 
  from step 2 is that this must use a different routine to copy
  the local stencil into the contiguous buffer, since these copies must 
  appropriately handle data from the  buffers.

  * Finally, compute the formula (<ref>) over the entire
  local subdomain, i.e., for each location (x_i,y_j,z_k) over the
  ranges i=0,…,n_xloc-1, j=0,…,n_yloc-1, and
  k=0,…,n_zloc-1.

We note that due to the high arithmetic intensity of the FD-WENO approach, 
each MPI task need not own a very large computational subdomain 
for the point-to-point communication to be completely overlaid by 
computation in step 2.




 §.§ Temporal discretization



After spatial semi-discretization, the original PDE system
(<ref>) may be written as an ODE initial-value problem,

    _t = (t,) + (t,),

where  contains the spatial semi-discretization of the solution vector 
, (t,) contains the spatial semi-discretization of the
terms ((,t)-∇·()) and (t,)
contains the spatial semi-discretization of the term ().
Here, we use the “S” superscript to denote the “slow” dynamical processes
(advection and externally-applied forces), and “F” to denote the
“fast” reaction processes.

For non-reactive flows in which =0, the initial value
problem is nonstiff, and is therefore solved using a
temporally-adaptive explicit Runge–Kutta method from ARKode's
 module.  While that use case is supported by our
demonstration code, we do not examine this “single physics” use case
here.

For problems involving chemical reactions, the multiphysics
initial-value problem (<ref>) typically exhibits multiple
temporal scales.  We therefore solve these problems using ARKode's
 module, that employs a third-order accurate
multirate infinitesimal step method for problems characterized by
two time scales <cit.>.  

Here, a single time
step to evolve (t_n-1) →(t_n-1+h^S), denoted by
_n-1→_n, for the full initial-value problem
(<ref>), proceeds according to the following algorithm:



  * set _1 = _n-1,

  * for i = 2,…,s+1:


  * define the “fast” initial condition: (t^S_n,i-1) = _i-1,

  * compute the forcing term
  
    = 1/c^S_i - c^S_i-1∑_j=1^i-1
        (A^S_i,j - A^S_i-1,j) (t^S_n,j, _j),


  * for τ∈ (t^S_n,i-1, t^S_n,i], solve the “fast”
  initial-value problem
  
    (τ) = (τ, ) + ,


  * set the new “slow” stage: _i = (t^S_n,i),


  * set the time-evolved solution: _n = _s+1,

where (A^S, b^S, c^S) correspond to the coefficients for an
explicit, s-stage “slow” Runge–Kutta method, A^S is
padded with a final row, A^S_s+1,j=b^S_j, and where 
t^S_n,j = t_n-1 + c^S_j h^S for j=1,…,s correspond to 
the “slow” stage times.  In this demonstration application, we use 
the default “KW3”  slow Runge–Kutta coefficients 
<cit.>.  For evolution of the “fast” problems 
(<ref>) above, we use a temporally-adaptive, 
diagonally-implicit Runge–Kutta method from ARKode's 
 module, namely the “ARK437L2SA” DIRK
method from <cit.>.

As we will describe in Section <ref> below when 
discussing our physical test problem, chemically-reactive flows exhibit 
fast transient behavior, so their stable evolution heavily depends on the
inherent robustness of temporally-adaptive integration.  However, the
primary purpose of this report is to document parallel performance of
the  and other algebraic solver enhancements
that have been recently added to SUNDIALS.  As such, we employ a
hybrid adaptive + fixed-step integration approach for the fast
time-scale subproblems (<ref>).  Specifically, we
partition the overall temporal domain (t_0,t_f] into two parts,
(t_0,t_0+h_t] and (t_0+h_t,t_f].  The first portion is considered
the “transient” time period, where chemical species exhibit
very fast dynamical changes, as they rapidly adjust from their
initial conditions to their slower (but still fast) solution
trajectories.  This time period is therefore evolved with
's temporal adaptivity enabled; the adaptivity
parameters employed for this phase of the simulations are provided in
Table <ref>.


The second (and typically much longer) time interval, (t_0+h_t,t_f],
is evolved using a fixed “fast” time step size, h^F = h^S/1000,
thereby bypassing 's built-in temporal adaptivity
approaches to produce a more predictable amount of work as the problem
is pushed to larger scales.  The challenge with using fixed time step 
sizes in such an application is that in fixed-step mode, any algebraic 
solver convergence failure becomes fatal, and causes the
entire simulation to halt.  Thus, the value of h_t must be chosen
appropriately, to balance the need for adaptivity-based robustness
during initial transient chemical evolution (i.e., larger h_t)
against the desire for a fixed amount of computational work
per MPI task when performing weak scaling studies (i.e., smaller
h_t).  We present the values used in the current studies in Section 
<ref>, when discussing the particular test problem 
used here.

Both evolution periods, (t_0,t_0+h_t] and (t_0+h_t,t_f], are
evolved using single calls to , although
the second period can be broken into a sequence of separate subperiods
so that solution statistics can be displayed,  and/or solution
checkpoint files can be written to disk.  However, since this study
focuses on overall solver performance, all such diagnostic and
solution output is disabled.





 §.§ Algebraic solvers



Since the slow time scale is currently treated explicitly, the only
algebraic solvers present in these calculations occur when evolving
each “fast” subproblem (<ref>).  As these subproblems
are evolved using diagonally implicit Runge–Kutta (DIRK) methods,
each stage may require the solution of a nonlinear algebraic system of 
equations.  We briefly outline the structure of these DIRK methods,
and then discuss the algebraic solvers used for these problems.

Considering the fast IVP

    (τ) = (τ, ) + ≡(τ,),

a σ-stage DIRK method with coefficients (A,b,c) evolves one
time step (τ_m-1) →(τ_m-1+θ), denoted for
short by _m-1→_m, via the algorithm: 


  * for i = 1,…,σ, solve for the “stages” _i that satisfy the equations
  
    _i = _m-1 + θ∑_j=1^i A_i,j(τ_m,j, _j),

  where τ_m,j = τ_m-1 + c_j θ,

  * compute the time-evolved solution
  
    _m = _m-1 + θ∑_i=1^σ b_i (τ_m,i, _i),


  * (optional) compute the embedded solution
  
    _m = _m-1 + θ∑_i=1^σb̃_i (τ_m,i, _i).



The solution of at most σ nonlinear algebraic systems (<ref>)
is required to compute each stage _i. Writing these systems in standard
root-finding form, for each fast substep we must solve σ
separate nonlinear algebraic systems:

    0 = (_i)
      ≡[_i - θ A_i,i(τ_m,i,_i)] - [_m-1 + θ∑_j=1^i-1 A_i,j(τ_m,j,_j)],
             i=1,…,σ.

where the first bracketed term contains the implicit portions of the
nonlinear residual, and the second bracketed term contains known
data.  We note that each _i ∈^n_x n_y n_z (5+n_c) can be
a very large vector.  However since = (τ, ) +, and (τ,) is just the spatially semi-discretized
version of the reaction function (), we see that although
equation (<ref>) is nonlinear, it only involves couplings
between unknowns that are co-located at each spatial location, 
(x_i,y_j,z_k).  Thus we may instead consider equation 
(<ref>) to be equivalent to a system of n_x n_y n_z
separate nonlinear systems of equations, each coupling only 
5+n_c unknowns.

We may therefore leverage this structure in a myriad of ways to
improve parallel performance.  At one extreme, we may break
apart the large nonlinear system (<ref>) into n_x n_y
n_z separate nonlinear systems, performing an independent Newton
iteration separately at each spatial location.  At a coarser 
level, we could instead break apart equation (<ref>) into
n_px n_py n_pz separate nonlinear systems, one per MPI task, 
so that each Newton iteration may proceed without parallel 
communication.  The relative merits of these choices (as well as 
intermediate options, e.g., breaking this into two-dimensional 
“slabs” or one-dimensional “pencils” of spatial cells) is not 
investigated here, but may be pursued in future efforts.  In this 
work, we utilize an even  coarser approach that solves the full 
nonlinear system of equations (<ref>)
using a single modified Newton iteration (the 
default), that couples all unknowns across the entire parallel
machine.  However, we exploit the problem structure by providing a
custom linear solver module that solves each MPI task-local
linear system independently.  Due to the structure of , the
Jacobian matrix J = ∂()/∂
is block-diagonal,

    J = [   J_1                  ;         J_2            ;                 ⋱      ;                   J_n_p ],

where each MPI task-local block J_p
∈^n_xlocn_ylocn_zloc(5+n_c), is itself block-diagonal,

    J_p = [                J_p,1,1,1                                                                           ;                                         J_p,2,1,1                                                  ;                                                                          ⋱                         ;                                                                            J_p,n_xloc,n_yloc,n_zloc ],

with each cell-local block J_p,i,j,k∈^5+n_c.  We leverage
the extreme sparsity of these MPI task-local Jacobian matrices J_p by
constructing each J_p matrix in compressed-sparse-row (CSR) format, and storing
them using standard  objects.  These are thus
automatically converted to Newton system matrices A = I-γ J_p
within 's generic direct linear solver interface.  We
then leverage the block-diagonal structure and solve the overall
Newton linear systems Ax=b using the  linear
solver module separately on each MPI task.

As this is an optimally-parallel direct linear solver for the
block-diagonal Newton linear systems, and since  itself
involves no parallel communication, the only “wasted” MPI
communication that occurs in our solution of each nonlinear algebraic
system (<ref>) is the computation of residual norms
()_WRMS to determine completion of the Newton
iteration.  In future studies we plan to supply a custom nonlinear
solver that will remove this extraneous communication.  We also plan
to explore alternate strategies that will further leverage the
block-diagonal structure (slabs, pencils, etc.), particularly once we
transition simulation of the chemical kinetics to GPU accelerators.





§ TEST DESCRIPTION



In this work, we consider the test problem of an advecting and
reacting primordial gas.  In addition to the five “standard” fluid
variables described in Section <ref>, we evolve 10
chemical species (i.e., n_c=10), for a total of 15 variables per
spatial location.  These species model the chemical behavior of a
low density primordial gas, present in models of the early universe
<cit.>.  This model consists of the species:


  *  – neutral atomic Hydrogen density

  * ^+ – positively-ionized atomic Hydrogen density

  * ^- – negatively-ionized atomic Hydrogen density

  * _2 – neutral molecular Hydrogen density

  * _2^+ – ionized molecular Hydrogen density

  *  – neutral atomic Helium density

  * ^+ – partially-ionized atomic Helium density

  * ^++ – fully-ionized atomic Helium density

  * ^- – free electron density

  * e_g – internal gas energy (proportional to temperature)

The full set of chemical rate equations that encode these reactions
comprise the reaction term, (), used in our code.  Both the
routine to evaluate this “right-hand side” function, as well as a
corresponding routine to evaluate its Jacobian in CSR format, are
provided by the Dengo software package <cit.>, a
source-code generation utility that translates from astropysical
chemical rate equations to  (or ) code that,
among other things, implements these routines and generates the
corresponding reaction rate lookup tables in HDF5 format
<cit.>.

We highlight the fact that the internal gas energy, e_g, may be
uniquely defined by the fluid fields, since

    e_t    = e_g + 1/2ρ(m_x^2 + m_y^2 + m_z^2)
    ⇔     
    
      e_g    = e_t - 1/2ρ(m_x^2 + m_y^2 + m_z^2),

and thus this reaction network only adds 9 new independent fields to
the simulation.  However, due to the multirate solver structure
described in Section <ref>, wherein
control over the simulation cleanly shifts between “slow” and
“fast” phases, we use both e_g and e_t to store the
“current” gas/total energy value.  Furthermore, by storing two
versions of the energy we may use our preferred 
data structure layout – five  objects for
the fluid fields, plus one MPI task-local  for the
chemistry fields.  We further note that this structure will enable
follow-on efforts in which the entire chemical network (data and
computation) are moved to GPU accelerators, through merely swapping
out our  object and enabling Dengo's
generation of GPU-enabled source code.

We initialize these simulations using a “clumpy” density field.
Here, the overall density field at any point ∈Ω is given
by the formula

    ρ(,t_0) = ρ_0(1 + 5 e^-20-_c^2 +
        ∑_i=1^10  n_p s_i e^-2(-_i/r_i)^2),

where ρ_0=1.67× 10^-22 g / cm^3 

is the background density, 
n_p is the total number of MPI tasks in the simulation (i.e., the number of “clumps” is proportional to the number of MPI tasks),

and
_c = (x_l+x_r/2, y_l+y_r/2, z_l+z_r/2) 
is the location of the center of the computational domain.  We choose
the remaining parameters from uniform random distributions in the
following intervals: 


  * _i ∈Ω, i.e, each clump is centered randomly within
  the domain,

  * r_i ∈ [3Δ x, 6Δ x] is the clump “radius,”
  i.e., each clump extends anywhere from 3 to 6 grid cells away from 
  its center, and

  * s_i ∈ [0, 5] is the clump “size,” i.e., each clump
  has density up to 5 times larger than the background density.


The simulation begins with a near-uniform temperature field, with only
a single higher-temperature region located in the clump at the domain
center:

    T(,t_0) = T_0(1 + 5 e^-20-_c^2),

where the background temperature is chosen to be T_0=10 K.

The chemical fields are initialized to values proportional to the
overall density, with:

    _2(,t_0)    = 10^-12ρ_0(,t_0)
    _2^+(,t_0)    = 10^-40ρ_0(,t_0)
    ^+(,t_0)    = 10^-40ρ_0(,t_0)
    ^-(,t_0)    = 10^-40ρ_0(,t_0)
    ^+(,t_0)    = 10^-40ρ_0(,t_0)
    ^++(,t_0)    = 10^-40ρ_0(,t_0)
    (,t_0)    = 0.24ρ_0(,t_0) - ^+(,t_0) - ^++(,t_0)
    (,t_0)    = ρ_0 - _2(,t_0) - _2^+(,t_0) - ^+(,t_0)
         - ^-(,t_0) - ^+(,t_0) - ^++(,t_0) - (,t_0)
    ^-(,t_0)    = ^+(,t_0)/1.00794 + ^+(,t_0)/4.002602 + 2^++(,t_0)/4.002602 - ^-(,t_0)/1.00794 + _2^+(,t_0)/2.01588
    
      e_g(,t_0)    = k_b T(,t_0) N(,t_0)/ρ(,t_0) (γ-1),

where k_b=1.3806488× 10^-16 g·cm^2/s^2/K is Boltzmann's constant, 

γ=5/3
is the ratio of specific heats, and the gas number density N(,t_0) is given by

    N(,t_0) = 5.988× 10^23(
         _2(,t_0)/2.01588 + _2^+(,t_0)/2.01588
        + ^+(,t_0)/1.00794 + ^-(,t_0)/1.00794
        + ^+(,t_0)/4.002602 + ^++(,t_0)/4.002602
        + (,t_0)/4.002602 + (,t_0)/1.00794).

Finally, we assume that the gas is initially static (i.e.,
𝐦(,t_0) = 0), and that 
no external forces are applied (i.e., (,t) = 0).

We perform the simulations on the spatial domain Ω =
[0,3.0857×10^30 cm]^3, and enforce reflecting boundary conditions
over ∂Ω.  Our base grid simulations (discussed below)
are computed over the temporal domain [0,10^11 s].   


Since the Dengo-supplied chemical network expects values in CGS units,
but our FD-WENO solver prefers dimensionless quantities, we
non-dimensionalize by using the scaling factors: 


  *  = 3× 10^70 g,

  *  = 3.0857× 10^30 cm, and

  *  = 10^11 s,

that correspond to the variable scaling factors


  *  = 1.0211× 10^-21 g/cm^3,

  *  = 3.1507× 10^-2 g/cm^2/s, and

  *  = 9.7223× 10^17 g/cm/s^2,

which we use to convert between “dimensional” and “dimensionless”
values as control is passed back-and-forth between fluid and chemistry
solvers.  Furthermore, we note that these choices of 
and  result in the normalized space-time domain
[0,1]^3×[0,1] – all discussion of time step sizes or the temporal 
domain in the remainder of this report refer to these quantities in 
dimensionless units.

As is typical for large-scale explicit simulations, we perform weak
scaling studies by examining the time required per slow step.
Thus as we refine the spatial discretization (e.g., Δ x →Δ x/2) and thus reduce the slow step size (e.g., h^S
→h^S/2), we shorten the overall time interval for each
simulation (e.g., [0,1]→[0,1/2]), thereby
guaranteeing a steady number of “slow” time steps.  However, since
the fast time scale calculations are performed implicitly and thus
have no resolution-based CFL stability restriction, we may choose
between two alternate options: 


  (a) maintain an essentially-constant fast step size (i.e., h^F
  → h^F), as would likely occur for chemical accuracy
  considerations alone, 

  (b) refine the fast step size in proportion to the slow step
  size (i.e., h^F →h^F/2), resulting in an essentially
  constant amount of work per MPI task per slow time step as the mesh
  is refined. 

Since the purpose of this report is to examine the parallel
scalability of the  and other general SUNDIALS
enhancements, we chose option (b) above, since option (a) would result 
in a decreasing number of fast time steps per simulation.  Thus the 
ensuing scalability results may be seen as a “worst-case” 
performance metric for scalability of multirate methods applied to 
physical problems of similar structure.   

We perform a standard explicit method weak scaling study for this
problem, wherein we increase the total number of nodes and total
problem size proportionately.  For each spatial mesh, we decrease the
“slow” step size h^S to maintain a constant CFL stability
factor (i.e., h^S ∝min(Δ x, Δ y, Δ z)).  In
each simulation, the initial “transient” phase runs over the time
interval (0,0.1], and the “fixed-step” phase runs for the remainder 
of the time interval.  We set a value h^F to
maintain a constant timescale separation ratio h^S/h^F = 1000; for the
transient phase h^F this is the maximum step size for the temporal
adaptivity approach; h^F then becomes the fixed step size for the 
latter portion of the calculation.  As stated above, we maintain an
essentially-constant computational effort per MPI task by shortening
the overall simulation time proportionately with h^S.  We note that
this causes the relative fraction of “fixed-step” versus
“transient” portions of the simulation to decrease at larger scales.
However, since we use the same h^F value for the fixed-step phase 
of the calculation as we use for the maximum allowed transient 
time step size, then as the spatial mesh is refined we find  
that an increasing fraction of transient chemical time steps use this 
maximum value, thereby providing a robust solution approach that in the 
limit achieves a fixed amount of work per MPI task as the mesh is refined.  
In Table <ref> we provide a
summary of these parameters for each problem size tested.


For each spatial mesh size, we performed two simulations.  The first
used two “newer” features added to the SUNDIALS 
API: 


  * Fused vector operations – for back-to-back
   operations that are performed in SUNDIALS'
  various solvers, we created a set of new  kernels
  that perform this multiple combination of operations per memory
  access, thereby increasing arithmetic intensity and reducing the
  number of function calls (or GPU kernel launches). 

  * Local reduction operations – since an
   is merely a vector that is comprised of a
  subset of other  objects, any operation requiring
  an MPI reduction (e.g., dot-product or norm) would naively
  result in multiple separate  calls.  
  We therefore created a set of new 
  routines that only perform the MPI task-local portion of these
  operations, waiting to call  until the
  overall accumulated quantity is available, and thereby reducing MPI
  latency effects. 

The second simulation performed at each spatial mesh size was run with
these newer features disabled, so that we could assess the benefits of
these newer features “at scale.” 





§ PERFORMANCE RESULTS



For each simulation we manually instrumented the code with 
profiling timers based on  for a variety of 
logical subsets of the code: 


  (a) Setup – this includes construction of ,
  , and  data structures, as well as
  construction of initial condition values.  We note that due to the
  sum over 10 n_p clumps in creation of the initial density field
  (<ref>), this component should exhibit slowdown as the
  problem size is increased.  We therefore do not include this 
  contribution when assessing the overall parallel efficiency. 

  (b) I/O – this includes all time spent in output of diagnostic
  information to  and  as the simulation
  proceeds, as well as all time spent in reading and writing HDF5
  files.  For the results that follow, most of these capabilities were
  disabled so that disk I/O would not affect our performance measurements. 

  (c) MPI – this includes all time spent in parallel
  communication by this demonstration application code.  This 
  does not include the MPI time spent in SUNDIALS'
   operations (e.g., dot-products and norms). 

  (d) Packing – this includes all time spent in packing the
  contiguous memory buffers (see steps 2 and 4 from Section
  <ref>) that are used in the FD-WENO flux
  reconstruction. 

  (e) FD-WENO – this includes all time spent in performing the
  FD-WENO flux reconstruction. 

  (f) Euler – this includes the time spent in all steps from
  Section <ref>, i.e., it should include the
  sum of (c)-(e) above, as well as the time required to evaluate the
  flux divergence (<ref>). 

  (g)  – this includes the time spent in (f)
  above, as well as all translation when converting between
  dimensional and dimensionless units for the fluid and chemical
  fields. 

  (h)  – this includes only the time spent in
  evaluating the ODE right-hand side for the Dengo-supplied chemical
  network. 

  (i)  – this includes only the time spent in
  evaulating the Jacobian of , and storing those values
  in CSR matrix format. 

  (j) Linear solver setup – this includes all time spent by our
  custom  implementation to factor its 
  block-diagonal linear system matrices; this is essentially just the 
  time spent in MPI task-local calls to the  
  “Setup” routine. 

  (k) Linear solver solve – this includes all time spent by our 
  custom  implementation to solve its block-diagonal 
  linear systems; it is essentially just the time spent in MPI task-local 
  calls to the  “Solve” routine. 

  (l) Overall “transient” simulation time – this includes all
  time spent in evolving the solution over the initial transient time
  interval, (t_0,t_0+h_t].  We note that this does not include the
  “Setup” time (a) above. 

  (m) Overall “fixed-step” simulation time – this includes all
  time spent in evolving the solution over the second fixed-fast-step
  time interval, (t_0+h_t,t_f]. 

  (n) Total simulation time – this includes both the
  “transient” and “fixed-step” simulation times, as well as the
  “Setup” time for the simulation, i.e., (a+l+m). 




We 
note that although these profilers only directly measure 
the time spent in the “multiphysics” portions of this simulation, we
may indirectly measure the overall amount of time spent in SUNDIALS
modules (, , , etc.)
by subtracting the amount of time spent in right-hand side, Jacobian,
and linear solver operations from the overall amount of time spent in
evolving the system, i.e.,

    SUNDIALS time = (l + m) - (g + h + i + j + k)

We further note that the vast majority of synchronization points
in this code occur in  reduction operations (e.g., 
dot-products and norms), and thus the “SUNDIALS time” includes all 
such synchronization points.  While some local synchronization occurs from 
nearest-neighbor communication for flux calculations, this is overlaid 
by flux computations over subdomain interiors, and is directly measured 
by the timers (c), (f), and (g) above. 

We provide weak-scaling results with these timers in Figure
<ref>.  While we collected 
data on a large number of code components, most of those timers required only 
trace amounts (less than 0.5%) of the total runtime: Setup, I/O, MPI, Packing, 
FD-WENO, Euler, , , and linear solver setup.  We 
have thus removed these measurements from this figure to 
focus on the time-intensive aspects of the 
code, as well as their parallel scalability.  Since each profiler resulted in 
slightly different times on each MPI task, for each curve we plot both the mean 
task time, as well as error-bars showing the minimum and maximum reported 
values.  However, since the `Sundials' times are only measured indirectly, 
these minimum and maximum values may be exaggerated since they incorporate the 
variability from other portions of the code (in addition to variability 
resulting from within SUNDIALS itself).  Thus the error bars for the `Sundials' curve in this Figure likely over-estimate the variability encountered by MPI tasks in the SUNDIALS infrastructure.



Additionally, in Table <ref> we provide the parallel efficiency of each simulation, in comparison with the 80-task “fused” simulation.



We first note that at the smallest scale tested (125× 100× 100 grid
with 80 MPI tasks), the fast chemical right-hand side routine accounted 
for almost 70% of the total runtime, followed by general SUNDIALS infrastructure 
(slightly over 20%), and linear system solves (slightly under 10%).  
Furthermore, of these the fast 
chemical right-hand side routine and linear solver times showed perfect 
weak scalability, while the time spent in SUNDIALS infrastructure increased 
by almost 42%, leading to run-time increase of about 11% for the overall simulation 
code (i.e., approximately 90% parallel efficiency compared to the smallest test).  
We additionally note that the variance 
in reported times increased at larger scales, particularly for the SUNDIALS
infrastructure, that showed up to 77% variance at 138,240 MPI tasks.  Interestingly,
the various MPI tasks showed essentially zero variability in the total simulation times 
reported.  This is likely due to the fact that these “total simulation” timers surrounded calls to 
; since the last thing this does before returning is compute an 
error norm (via ), it effectively synchronizes all MPI tasks just 
prior to stopping the timer.
Finally, we note that the “fused” and “unfused” versions showed no
statistically-significant differences.



§ CONCLUSIONS



While we expected the general strong performance results shown in Section
<ref>, some of these came as more of a surprise than others.

From an application viewpoint, we note that our approach for interleaving 
communication and computation for the advective fluxes proved sufficient, with 
all  directly-measured “MPI” timings remaining at essentially zero for all 
scales tested.  Moreover, due to the multirate structure of this application
problem and testing setup, the advective portion of the right-hand side was 
essentially “free” in comparison to the chemical kinetics that dominated 
the simulation time.

From a SUNDIALS-infrastructure viewpoint, we posit that the slowdown shown in
the SUNDIALS infrastructure was entirely due to synchronization-induced 
latency in  calls.  This conclusion is based on two 
observations.  First, while the “fused” version reduced the number 
of these calls by 80%, it left the frequency of these synchronization
points essentially unchanged, and thus the near-identical runtimes for both 
indicate that the additional  calls in the “unfused” 
version do not increase the overall global synchronization of the code.  
Moreover, the large variance reported 
for SUNDIALS infrastructure among MPI tasks indicates that some tasks spent 
considerably more time waiting at these synchronization points than others.
However, since we could only indirectly measure the time spent in SUNDIALS 
routines, this variance may be exaggerated.

Based on the above, we anticipate that this code will tremendously benefit 
from three planned investigations.  First, we plan to extend this implementation to
construct a custom  for the fast time scale.  As 
discussed in Section <ref>, this can be constructed 
to remove all  calls at the algebraic solver level, in 
turn reducing over 98% of the overall  calls from the 
“fast” time scale (for the 138,240 MPI task run, there were 61,195 Newton 
iterations and 1,007 fast time steps; 61195/(61195+1007) ≈ 0.984), thereby 
effectively removing nearly all global synchronization from these simulations.

Second, we plan to port the chemical network data, the fast right-hand side 
computations, the fast Jacobian construction routine, and the fast-time scale 
algebraic solvers to GPU accelerators.  



We note that even at the largest scales tested here
(1500 × 1200 × 1200 grid on 138,240 MPI tasks), evolution 
of these arithmetically-intense “fast” chemical processes required over 
70% of the total runtime (fast RHS + lsolve).  Moreover, 
the types of calculations performed in this module should be amenable to GPU
architectures, and thus their conversion could likely result in 
significant performance improvements over the results shown here.

Third, instead of requiring application codes to only indirectly measure
the time spent in SUNDIALS, and thus lumping all of SUNDIALS' various
actions (plus any unmeasured work) into a single performance metric, 
future performance studies with SUNDIALS would benefit tremendously from 
direct measurements of SUNDIALS performance, notably if these 
measurements were broken apart into logical units (e.g., vector reductions, 
vector arithmetic, integrator infrastructure, nonlinear solvers, linear 
solvers, preconditioners, etc.).  We are thus investigating inclusion of 
direct measurements of SUNDIALS performance through a tool such as Caliper
<cit.>.  Additionally, such tools could track more than just 
runtime, enabling measurement of system-level counters.



§ ACKNOWLEDGEMENTS

This research was supported by the Exascale Computing Project (ECP), Project Number: 17-SC-20-SC, a collaborative effort of two DOE organizations – the Office of Science and the National Nuclear Security Administration, responsible for the planning and preparation of a capable exascale ecosystem, including software, applications, hardware, advanced system engineering and early testbed platforms, to support the nation’s exascale computing imperative.

Support for this work was also provided through the Scientific Discovery through Advanced Computing (SciDAC) project “Frameworks, Algorithms and Scalable Technologies for Mathematics (FASTMath),” funded by the U.S. Department of Energy Office of Advanced Scientific Computing Research and National Nuclear Security Administration. 

This research used resources of the Oak Ridge Leadership Computing Facility, which
is a DOE Office of Science User Facility supported under Contract DE-AC05-
00OR22725.

This work work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344. Lawrence Livermore National Security, LLC.  LLNL-TR-791538.

This document was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes.

siam


