CS 5220

Applications of Parallel Computers

Distributed parameter models

Please click the play button below.

Types of PDEs

Type	Example	Time?	Space dependence?
Elliptic	electrostatics	steady	global
Hyperbolic	sound waves	yes	local
Parabolic	diffusion	yes	global

All right. If you're going to talk about PDEs at a cocktail party, you should know the names of the main classes of PDEs that everyone around you is talking about.

First, there are elliptic PDEs. These come up in equilibrium problems, like finding the electrostatic potential for a system. They're steady-state problem, so no real time dependence; but the solution at one point in space depends on the behavior everywhere else in space, at least within the domain. So we say that this class of problems has global dependence in space.

Second, there are hyperbolic PDEs. These are the equations that describe waves and signals and the like. The usual model problem people talk about is propagation of sound. This is a time-dependent problem, so it's more complicated then elliptic problems in that way. But waves travels through space at a finite speed, so points that are far from each other can be dealt with more-or-less independently, at least over short periods.

Yes, yes, sometimes we do frequency-domain analysis of sound. Who pointed that out? Oh, I guess it's just me here. Moving on, then...

The third type of PDE is parabolic equations. Think diffusion processes here, like the diffusion of heat. These problems are time-dependent, but the equations don't give us anything like a speed of sound. Everything affects everything else instantaneously, though maybe not by a lot.

The differences between these types of equations turn out to matter a lot when it comes time to design - and parallelize - numerical solution algorithms. We'll sketch why in the rest of this deck.

And now, you're ready for that cocktail party! Or you would be, if we were going to cocktail parties any more. And if you frequented the types of cocktail parties attended by numerical PDE folks.

OK, let's just move on with the lecture, shall we?

Types of PDEs

Different types involve different communication:

Global dependence \(\implies\) lots of communication
(or tiny steps)
Local dependence from finite wave speeds;
limits communication

Example: 1D heat equation

Consider flow (e.g. of heat) in a uniform rod

Heat (\(Q\)) \(\propto\) temperature (\(u\)) \(\times\) mass (\(\rho h\))
Heat flow \(\propto\) temperature gradient (Fourier’s law)

Our first example will be heat flow in a uniform rod (that is, a 1D medium where the heat capacity is constant).

In this setting, the amount of heat in a region is proportional to the temperature times the mass. Heat flows from hotter regions to colder regions, tending to even out over time. We sometimes introduce this in terms of Fourier's law, which says that the heat flow is proportional to the temperature gradient. But I learned another derivation when I first saw this stuff, in terms of Newton's law of cooling.

Newton's law of cooling says that we lose heat across an interface at a rate proportional to the temperature jump across the interface (and the surface area of that interface). This is the way we figure out how long a body has been dead, for example, where by "we" I mean somebody completely other than me - I've never been called on to do this in person.

So suppose we think of our long, uniform rod as a sequence of short, uniform rods, each with roughly constant temperature. The rate at which heat moves from one of these short segments to its neighbors is determined by Newton's law of cooling. If we let the length of these short segments go to zero, in the limit we recover Fourier's law.

Example: 1D heat equation

Consider flow (e.g. of heat) in a uniform rod

Heat (\(Q\)) \(\propto\) temperature (\(u\)) \(\times\) mass (\(\rho h\))
Heat flow \(\propto\) temperature gradient (Fourier’s law)

\[\begin{aligned} \frac{\partial Q}{\partial t} \propto h \frac{\partial u}{\partial t} &\approx C \left[ \left( \frac{u(x-h)-u(x)}{h} \right) + \left( \frac{u(x)-u(x+h)}{h} \right) \right] \\ \frac{\partial u}{\partial t} &\approx C \left[ \frac{u(x-h)-2u(x)+u(x+h)}{h^2} \right] \rightarrow C \frac{\partial^2 u}{\partial x^2} \end{aligned}\]

Spatial discretization

Heat equation with \(u(0) = u(1) = 0\) \[\frac{\partial u}{\partial t} = C \frac{\partial^2 u}{\partial x^2}\]

Spatial discretization

Spatial semi-discretization: \[\frac{\partial^2 u}{\partial x^2} \approx \frac{u(x-h)-2u(x)+u(x+h)}{h^2}\]

Spatial discretization

Yields a system of ODEs \[\frac{du}{dt} = C h^{-2} (-T) u = -C h^{-2} \begin{bmatrix} 2 & -1 & & & & \\ -1 & 2 & -1 & & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_{n-2} \\ u_{n-1} \end{bmatrix}\]

Explicit time stepping

Approximate PDE by ODE system (“method of lines”): \[\frac{du}{dt} = C h^{-2} T u\] Now need a time-stepping scheme for the ODE.

Explicit time stepping

Simplest scheme is Euler: \[u(t+\delta) \approx u(t) + u'(t) \delta = \left( I - C \frac{\delta}{h^2} T \right) u(t)\]
Taking a time step \(\equiv\) sparse matvec with \(\left( I - C \frac{\delta}{h^2} T \right)\)
This may not end well...

Explicit data dependence

Plot of information propagation for explicit heat equation time stepper

Nearest neighbor interactions per step \(\implies\)
finite rate of numerical information propagation

Explicit time stepping in parallel

for t = 1 to N
  communicate boundary data ("ghost cell")
  take time steps locally
end

Overlapping communication with computation

for t = 1 to N
  start boundary data sendrecv
  compute new interior values
  finish sendrecv
  compute new boundary values
end

Batching time steps

for t = 1 to N by B
  start boundary data sendrecv (B values)
  compute new interior values
  finish sendrecv (B values)
  compute new boundary values
end

The problem with the scheme we've discussed so far is that we can do a lot of work in the time it takes to send one message - so much work, in fact, that we might still be unable to hide the cost of the communication underneath of our computation. But a big part of the communication overhead comes not from the volume of data that we are sending - the bandwidth costs - but from message latency. And we can avoid paying that latency cost so often by batching up the communication events. Instead of sending boundary data before every time steps, we send boundary data every B time steps, then operate independently. If we're going to go for B steps, we need to send B layers of boundary data in that message. Nothing changes the fundamental fact that we're propagating one boundary layer per time step on average; we are just changing the frequency with which we communicate, and thus the number of latency times we have to hide.

Explicit pain

Plot of instability in explicit heat equation stepper

Unstable for \(\delta > O(h^2)\)!

All right, now that we're feeling very clever about our skills with organizing communication, let's look what happens when we apply this technique to the heat equation. As we foreshadowed, the result doesn't look good. This is a 1D problem, with time running from back to front and space running from left to right. And, as you can see, the nice smooth temperature profile that we started with develops oscillations that grow in magnitude as time proceeds. This is what I meant when I said the simulation blows up. This is not a physical phenomenon; it's just a numerical artifact that comes from taking time steps that are too long for the discretization. Unfortunately, this will happen in general for the heat equation unless we take time steps smaller than the square of the spatial step.

Implicit time stepping

Backward Euler uses backward difference for \(d/dt\) \[u(t+\delta) \approx u(t) + u'(t + \delta t) \delta\]
Time step \(\equiv\) apply \(\left( I + C \frac{\delta}{h^2} T \right)^{-1}\)
No time step restriction for stability (good!)
But each step involves linear solve (not so good!)
- Good if you like numerical linear algebra?

As we've mentioned before, we can avoid the instability that we just saw by switching from an explicit time-stepper to an implicit scheme, something like backward Euler. Where forward Euler involves multiplying by a tridiagonal matrix at every step, backward Euler involves solving a tridiagonal linear system of equations at every step. With backward Euler, we have no time step restrictions for stability, though we do still need to control the time step in order to maintain good accuracy. But each step involves solving a linear system, which is maybe a bit more complicated than doing a sparse matvec. So this is a complication, though it's good news for any numerical linear algebra fans following along from home.

Explicit and implicit

Explicit:

Propagates information at finite rate
Steps look like sparse matvec (in linear case)
Stable step determined by fastest time scale
Works fine for hyperbolic PDEs

So what have we learned? Explicit time-steppers like forward Euler tend to propagate information at a finite rate through a numerical mesh, at least when we are using a finite difference discretization like the one we've discussed so far. In the linear case, taking a time step is just a sparse matvec. But we get into trouble if the physics can move information around our domain faster than the numerical method can keep up. The stable step is limited by the fastest time scale, and that fastest time scale is effectively infinite for parabolic differential equations like the heat equation. On the other hand, this approach works fine for hyperbolic PDEs like the wave equation, where the physics has a natural speed of information propagation to it. We just need to make sure that our numerics keeps up with the physics by not letting the time steps get too long. If you hang around with numerical PDE people, you'll sometimes hear this referred to as the CFL condition.

Explicit and implicit

Implicit:

No need to resolve fastest time scales
Steps can be long... but expensive
- Linear/nonlinear solves at each step
- Often these solves involve sparse matvecs
Critical for parabolic PDEs

In contrast to explicit time steppers, implicit time steppers (like backward Euler) don't need to resolve the fastest time scales in the system. This means that we can generally take longer steps than we would with an explicit stepper, and not get into so much trouble with numerical instability (though we still have to think about the role of the time step in overall solution accuracy). But we may pay more per step than we do with an explicit method, because we need to solve a system of equations at every step. And our favorite solution algorithms involve... sparse matvecs again! You can't get away from them. Of course, our favorite solver algorithms also usually involve some other ingredients, but we'll get to that later in the course.

When we are solving parabolic PDEs like the heat equation, we really have very little choice. Implicit time steppers are the only viable option. For hyperbolic PDEs, we can go explicit or implicit, though explicit methods are often favored.

Poisson problems

Consider 2D Poisson \[-\nabla^2 u = -\frac{\partial^2 u}{\partial x^2} - \frac{\partial^2 u}{\partial y^2} = f\]

Prototypical elliptic problem (steady state)
Similar to a backward Euler step on heat equation

So far, we've been talking a lot about parabolic and hyperbolic equations. But at the start of the slide deck, we also mentioned elliptic equations, usually associated with steady-state problems. The prototypical example of an elliptic problem is the Poisson equation, shown here.

You might notice that the Poisson equation looks a lot like the type of equation we would see in a backward Euler step on the heat equation. Turns out that figuring out how to solve Poisson fast is also useful for figuring out how to efficiently implement implicit solvers for the heat equation.

We considered just one spatial dimension in the case of the heat equation, but had to deal with time. For Poisson, we don't have to worry about time, so let's consider two spatial dimensions.

Poisson problem discretization

\(u_{i,j} = h^{-2} \left( 4u_{i,j}-u_{i-1,j}-u_{i+1,j}-u_{i,j-1}-u_{i,j+1} \right)\)

\[L = \left[ \begin{array}{ccc|ccc|ccc} 4 & -1 & & -1 & & & & & \\ -1 & 4 & -1 & & -1 & & & & \\ & -1 & 4 & & & -1 & & & \\ \hline -1 & & & 4 & -1 & & -1 & & \\ & -1 & & -1 & 4 & -1 & & -1 & \\ & & -1 & & -1 & 4 & & & -1 \\ \hline & & & -1 & & & 4 & -1 & \\ & & & & -1 & & -1 & 4 & -1 \\ & & & & & -1 & & -1 & 4 \end{array} \right]\]

We might think about the discrete solution as a collection of values on a 2D grid, but for setting up solvers, it's worthwhile thinking about a 1D indexing scheme. Let's consider what happens when we lay out the mesh points in column major order. In that case, we can write down the finite difference approximation in terms of a sparse matrix-vector product with a matrix like the one shown here. This is for a three-by-three mesh, because otherwise we would run out of room on the slide! The matrix is block tridiagonal, with tridiagonal blocks on the main diagonal. So it's a little more complicated than the tridiagonal matrix T we saw before, but related.

If this was my matrix computations class, this is the part of the lecture where I'd start going on and on about Kronecker products and the vec operator. But... it's not my matrix computations class. So let's move on for now.

Poisson solvers in 2D/3D

\(N = n^d =\) total unknowns

Ref: Demmel, Applied Numerical Linear Algebra, SIAM, 1997.

Remember: best MFlop/s \(\neq\) fastest solution!

People have thought a lot about numerical methods for solving the Poisson equation, or things that look like the Poisson equation. My presentation is stolen from Jim Demmel's book on numerical linear algebra (with full disclosure that Jim was one of my advisors in grad school).

I'm going to show you the asymptotic complexity results in the next couple slides. You might object that this goes against what I told you earlier about not putting too much faith in asymptotic complexity. But remember, too, that performance is about both the rate at which work gets done and the total amount of work. The simpler methods may get better flop rates than the more complicated methods, but the methods with better asymptotic complexity still win on big problems. And when the number of unknowns is the square (or cube in 3D) of the number of points per side of the mesh, our problems get big pretty fast.

Poisson solvers in 2D/3D

Method	Time	Space
Dense LU	\(N^3\)	\(N^2\)
Band LU	\(N^2\) (\(N^{7/3}\))	\(N^{3/2}\) (\(N^{5/3}\))
Jacobi	\(N^2\)	\(N\)
Explicit inv	\(N^2\)	\(N^2\)

We usually start with talking about pretty general methods that work with big classes of linear systems. Methods based on Gaussian elimination solve all types of linear systems, but they are expensive - N^3 space, where N is the number of mesh points in this case. You aren't going to use these for even toy problems in 3D; it just costs too much. These methods use a lot of space, too. We do a little better with band variants that use the fact that all the entries far from the diagonal are zero, but it's still expensive.

At the same complexity level as band Gaussian elimination, there's the simplest iterative method we're going to discuss: Jacobi's method. That turns out to take N^2 time to reduce the error by a constant amount.

We've also written down computing an explicit inverse (or, if you are a PDE person, working with a solver based on Green's functions). Because we know the Green's function for this type of problem, we can do this in N^2 time. It gets more complex when we can't write down the Green's function so easily.

You don't need to worry if you've never heard of a Green's function before, by the way. The point is that this is still a pretty expensive thing to do.

Poisson solvers in 2D/3D

Method	Time	Space
CG	\(N^{3/2}\)	\(N\)
Red-black SOR	\(N^{3/2}\)	\(N\)
Sparse LU	\(N^{3/2}\)	\(N \log N\) (\(N^{4/3}\))
FFT	\(N \log N\)	\(N\)
Multigrid	\(N\)	\(N\)

The faster methods for solving Poisson problems rely more and more on the structure of the PDE. The method of conjugate gradients and red-black SOR iteration are general-purpose iterative solvers, but the analysis of their rate of convergence on Poisson discretizations depends a lot on the nature of the Poisson problem. We can get these to run in O(N^3/2) time on this problem, and we'll talk about that later. We can also solve 2D Poisson in O(N^3/2) using a sparse LU factorization; there, we don't care so much that we are solving Poisson, but we do care that the sparsity pattern is associated with a nearest-neighbor mesh in 2D. Finally, we can get log-linear time using Fourier transforms, or even linear time using multigrid methods. But those are really tied to this type of linear system.

We will talk about many of these types of linear solvers later in the semester.

General implicit picture

Implicit solves or steady state \(\implies\) solving systems
Nonlinear solvers generally linearize
Linear solvers can be
- Direct (hard to scale)
- Iterative (often problem-specific)
Iterative solves boil down to matvec!

All right. So if we are doing implicit time stepping or solving steady-state PDEs like Poisson, we end up solving large linear systems of equations. And if we're doing more complicated nonlinear problems, we solve nonlinear systems of equations; but the nonlinear solvers we use are things like Newton's iteration, which involves a linear solve at every step. So for a lot of reasons, we are going to care about linear solvers.

As I alluded to a moment ago, there are two general classes of linear solvers. Direct methods like Gaussian elimination just directly solve the problem. They are great when they apply, but they cost too much to be used for lots of big PDE solvers. So for solving PDEs, we often focus on iterative methods, which take an initial guess and incrementally make it better and better until we decide that the level of error is acceptable. Making those iterative methods converge fast often involves some problem-specific thinking. But there's also a piece that is crucial to the step-by-step performance that we can think about in a somewhat problem-independent way: and that piece is the sparse matvec operation.

PDE solver summary

Can be implicit or explicit (as with ODEs)

Explicit (sparse matvec) — fast, but short steps?
- works fine for hyperbolic PDEs
Implicit (sparse solve)
- Direct solvers are hard!
- Sparse solvers turn into matvec again

OK, let's wrap up our lightning intro to numerical PDEs. For time-dependent PDEs, we can use either implicit or explicit solvers, just like with ODEs. We often choose the method based on the type of PDE: implicit for parabolic equations, explicit for hyperbolic (though we could also do implicit methods there, if we wanted).

Implicit methods generally require solving sparse linear or nonlinear systems of equations. Direct solution algorithms like Gaussian elimination don't necessarily scale well to problems the size of PDE meshes, so we are likely to turn to iterative solvers. And both iterative solvers and explicit time steppers lean pretty heavily on sparse matvecs, at least for linear PDEs.

PDE solver summary

Differential operators turn into local mesh stencils

Matrix connectivity looks like mesh connectivity
Can partition into subdomains that communicate only through boundary data
More on graph partitioning later

PDE solver summary

Not all nearest neighbor ops are equally efficient!

Depends on mesh structure
Also depends on flops/point

CS 5220

Applications of Parallel Computers

Distributed parameter models

Types of PDEs

Types of PDEs

Example: 1D heat equation

Example: 1D heat equation

Spatial discretization

Spatial discretization

Spatial discretization

Explicit time stepping

Explicit time stepping

Explicit data dependence

Explicit time stepping in parallel

Overlapping communication with computation

Batching time steps

Explicit pain

Implicit time stepping

Explicit and implicit

Explicit and implicit

Poisson problems

Poisson problem discretization

Poisson solvers in 2D/3D

Poisson solvers in 2D/3D

Poisson solvers in 2D/3D

General implicit picture

PDE solver summary

PDE solver summary

PDE solver summary

Onward