CS 5220

Applications of Parallel Computers

Lumped parameter models

Please click the play button below.

Lumped parameter simulations

Examples include:

SPICE-level circuit simulation
- nodal voltages vs. voltage distributions
Structural simulation
- beam end displacements vs. continuum field
Chemical concentrations in stirred tank reactor
- concentrations in tank vs. spatially varying concentrations

What do we mean when we talk about lumped parameter simulations? Basically, we mean systems of differential equations. These types of systems occur all over engineering: think circuit simulations, structural simulations, or models of chemical reactors. We call these "lumped" models because we make simplifying assumptions about physical quantities like voltages, displacement fields, or chemical concentrations. With these assumptions, we can describe the fields in terms of relatively few variables. In contrast, distributed-parameter models work with more faithful representations of the continuous fields - though everything is still on the computer, so we are only ever really going to deal with a finite number of variables when we get down into the code.

The terms "lumped parameter" and "distributed parameter" mostly come from certain corners of engineering; if you prefer, I won't mind if you call these ODE and PDE models.

Lumped parameter simulations

Typically ordinary differential equations (ODEs)
Constraints: differential-algebraic equations (DAEs)

Often (not always) sparse.

Sparsity

Consider ODEs $x' = f(x)$ (special case: $f(x) = Ax$)

Dependency graph: edge $(i,j)$ if $f_j$ depends on $x_i$
Sparsity means each $f_j$ depends on only a few $x_i$
Often arises from physical or logical locality
Corresponds to $A$ being sparse (mostly zeros)

When we talk about sparsity in this setting, we are really talking about a dependency graph between variables. In this graph, each variable is represented by a node, and we have an edge between i and j if the derivative of x_j depends on x_i. In the special case of a linear system of ODEs, sparsity of the system refers to sparsity of the system matrix A (or the property that most of the entries of A are zero). But we can talk about sparsity of nonlinear equations, too.

Why do we expect to see sparsity in simulations of physics or of engineered systems? Because of locality! In the case of physics, we expect things that are close together to interact more than things that are far away. In the case of engineered systems, we often try to build the system so that various things interact in rather controlled ways. Either way, sparsity is pretty common across many systems we'll deal with.

Sparsity and partitioning

Sparsity plot and associated graph, partitioned

Want to partition sparse graphs so that

Subgraphs are same size (load balance)
Cut size is minimal (minimize communication)

We’ll talk more about this later.

Types of analysis: Static

Consider $x' = f(x)$ (special case: $f(x) = Ax + b$).

Might want: $f(x_*) = 0$

Boils down to $Ax = b$ (e.g. for Newton-like steps)
Can solve directly or iteratively
Sparsity matters a lot!

Types of analysis: Dynamic

Consider $x' = f(x)$ (special case: $f(x) = Ax + b$).

Might want: $x(t)$ for many $t$

Involves time stepping (explicit or implicit)
Implicit methods involve linear/nonlinear solves
Need to understand stiffness and stability issues

In addition to finding steady-state solutions, we often want to run simulations that see what happens to the system from a given initial state. Because computers are not so good at dealing with continuous variables like time, these dynamic simulations approximate the problem by taking discrete time steps. Depending on the nature of the problem, we might use an explicit method or an implicit method; the latter involves solving a linear or nonlinear system of equations at every step. We choose between implicit and explicit methods based on the so-called "stiffness" of the problem: for stiff problems, explicit time-steppers blow up unless we take very small time steps. Differential-algebraic equations are infinitely stiff, so we always use implicit methods for time-stepping them.

Types of analysis: Modal

Consider $x' = f(x)$ (special case: $f(x) = Ax + b$).

Might want: eigenvalues of $A$ or $f'(x_*)$

Explicit time stepping

Example: forward Euler $$ x_{k+1} = x_k + (\Delta t) f(x_k) $$
Next step depends only on earlier steps
Simple algorithms
May have stability/stiffness issues

When we talk about numerically solving differential equations in introductory courses in differential equations, we usually start by describing Euler's method. This is based on a Taylor expansion of the function f around the state x_k at time t_k. In pictures, we draw a vector field from each point, and approximate the solution by following the vector field at every step. The resulting algorithm is simple and explicit: we don't need to solve an equation to get x_{k+1}, we only need to evaluate f at the previous position x_k.

Unfortunately, for many interesting systems, forward Euler is unstable unless we take rather short steps. The same is also true of more sophisticated explicit time-stepping schemes.

Implicit time stepping

Example: backward Euler $$ x_{k+1} = x_k + (\Delta t) f(x_{k+1}) $$
Next step depends on itself and on earlier steps
Algorithms involve solves — complication, communication!
Larger time steps, each step costs more

A common kernel

In all cases, lots of time in sparse matvec:

Iterative linear solvers: repeated sparse matvec
Iterative eigensolvers: repeated sparse matvec
Explicit time marching: matvecs at each step
Implicit time marching: iterative solves (involving matvecs)

We need to figure out how to make matvec fast!

Somewhat surprisingly - or not, if you're a linear algebra fan - sparse linear algebra algorithms sit at the heart of almost all the types of analysis methods we've described. More specifically, the fundamental building block for a lot of algorithms with large sparse matrices is multiplication of a sparse matrix by a vector (also called sparse matvec). This building block appears when we're solving linear systems, computing eigenvalues and eigenvectors, time-stepping with explicit methods, or solving systems of equations in a static analysis or an implicit time-stepper.

All of this is to say: we really care about the performance of sparse matvec, and will spend serious time talking about it in later parts of the class.

An aside on sparse matrix storage

Sparse matrix $\implies$ mostly zero entries
- Can also have “data sparseness” — representation with less than $O(n^2)$ storage, even if most entries nonzero
Could be implicit (e.g. directional differencing)
Sometimes explicit representation is useful
Easy to get lots of indirect indexing!
Compressed sparse storage schemes help

Example: Compressed sparse row

Illustration of compressed sparse row format

This can be even more compact:

Could organize by blocks (block CSR)
Could compress column index data (16-bit vs 64-bit)
Various other optimizations — see OSKI

One example of a compressed storage format for sparse matrices is the so-called compressed sparse row (or CSR) format. We've drawn a picture here to illustrate how it works. The format consists of three arrays: parallel data and column index arrays that store the value and column of each nonzero in the matrix, listed in row-major order; and a pointer array that indicates the offset of each row in the data and column arrays. This turns out to be a pretty efficient general-purpose format for sparse matrices where we want to compute fast matvecs.

One can try to be fancier, of course, particularly if there is some internal structure. We might organize the nonzero entries into blocks, or compress the column index array. These types of tweaks, along with several others, are the key building blocks for codes like OSKI, which produces tuned sparse matrix-vector products for specific matrices.

Summary

ODE and DAE models widely used in engineering
Different analyses: static, dynamic, modal
Sparse linear algebra is often key