1. No category

Introduction to Numerical Methods for Solving
Partial Differential Equations
Benson Muite
benson.muite@ut.ee
http://kodu.ut.ee/˜benson
http://en.wikibooks.org/wiki/Parallel_Spectral_Numerical_Methods
25 March 2015
Outline
• Review
• Aim
• The heat equation
• Methods for solving the heat equation
• Fourier Series
• The Allen Cahn Equation
• Summary
• References
Ordinary Differential Equations: Numerical Schemes
• Forward Euler method
y n+1 − y n
= f yn
δt
• Backward Euler method
y n+1 − y n
= f y n+1
δt
• Implicit Midpoint rule
n+1
y
+ yn
y n+1 − y n
=f
δt
2
• Crank Nicolson Method
f y n+1 + f (y n )
y n+1 − y n
=
δt
2
• Other Methods: Runge Kutta, Adams Bashforth, Backward
differentiation, splitting
The Heat Equation
• Used to model diffusion of heat, species,
• 1D
• 2D
• 3D
∂u
∂2u
=
∂t
∂x 2
∂u
∂2u ∂2u
=
+
∂t
∂x 2 ∂y 2
∂u
∂2u ∂2u ∂2u
=
+
+
∂t
∂x 2 ∂y 2
∂z 2
• Not always a good model, since it has infinite speed of
propagation
• Strong coupling of all points in domain make it
computationally intensive to solve in parallel
¨
Linear Schrodinger
Equation
• Used to model quantum mechanical phenomena and often
appears in simplified wave propagation models
• 1D
i
• 2D
i
• 3D
i
∂u
∂2u
=
∂t
∂x 2
∂2u ∂2u
∂u
=
+
∂t
∂x 2 ∂y 2
∂u
∂2u ∂2u ∂2u
=
+
+
∂t
∂x 2 ∂y 2
∂z 2
• Looks like a heat equation with imaginary time
• Strong coupling of all points in domain make it
computationally intensive to solve in parallel
The Wave Equation
• Used to model propagation of sound, light
• 1D
∂2u
∂2u
=
∂t 2
∂x 2
• 2D
• 3D
∂2u ∂2u
∂2u
=
+
∂t 2
∂x 2 ∂y 2
∂2u
∂2u ∂2u ∂2u
=
+
+
∂t 2
∂x 2 ∂y 2
∂z 2
• Has finite speed of propagation
• Finite signal propagation speed sometimes useful in
parallelization since there is no coupling between grid
points that are far apart, hence smaller communication
requirements
• If propagation speed is very fast, communication
requirements still important
The Poisson Equation
• Used to model static deflection of a drumhead
• 1D
f (x) =
• 2D
f (x, y) =
• 3D
f (x, y, z) =
∂2u
∂x 2
∂2u ∂2u
+
∂x 2 ∂y 2
∂2u ∂2u ∂2u
+
+
∂x 2 ∂y 2
∂z 2
• No time dependence, but also arises in time discretizations
of time dependent partial differential equations
Burgers’ Equation
• Simple model for gas dynamics, also traffic
• Inviscid form (can develop shocks)
ut = −uux
• Viscous form
ut = −uux + uxx
Navier-Stokes Equations
• Models motion of fluids
• Consider incompressible case only (low speed motion
relative to speed of sound)
•
ρ
∂u
+ u · ∇u
= −∇p + µ∆u
∂t
∇ · u = 0.
• p pressure, µ viscosity, ρ, density
• 2D u(x, y ) = (u(x, y), v (x, y))
• 3D u(x, y , z) = (u(x, y, z), v (x, y, z), w(x, y, z))
Maxwells Equations
• Model propagation of electromagnetic waves
•
~t − ∇ × H
~ =0
D
•
~t + ∇ × E
~ =0
B
• with the relation between the electromagnetic components
given by the constitutive relations:
•
~ = εo (x, y, z)E
~
D
•
~ = µo (x, y, z)H
~
B
Summary on Models
• Many diverse models
• In 2D and 3D, parallel computing is very useful for getting
numerical solutions in reasonable time
• Typical applications are in physics, chemistry, biology and
engineering
• Related models also appear in social sciences, though
usefulness for predictability of real world data is less clear
Aim to simulate this,
•
https://www.flickr.com/photos/kitware/2293740417/in/pool-paraview/
• Visualization around Formula 1 Race Car by Renato N.
Elias, Rio de Janeiro, Brazil
or this
•
https://www.flickr.com/photos/kitware/2294528826/in/pool-paraview/
• Visualization around Formula 1 Race Car by Renato N.
Elias, Rio de Janeiro, Brazil
Maybe even this
• https://youtu.be/FTlJ9l9KTkw
• https://youtu.be/MsDgw8_cb90
• https://youtu.be/Su8qzC4HHVs
A simple model: The Heat Equation
The Heat Equation
• How does the temperature of the rod vary with time?
The Heat Equation: Model 1
• Middle of rod is initially hot due to previous heating (eg.
oven in house that has just been turned off)
• Temperature on the two sides is 0 (winter and cold outside
the house)
• Assume discrete uniformly space time, and discrete space
with molecules at each coordinate point
• Assume temperature in bar is directly proportional to heat,
with the same proportionality constant throughout the bar
• At next time step, a molecule will
• with probability 0.5 not transfer its current heat
• with probability 0.25 transfer half of its heat to its left
neighbor
• with probability 0.25 transfer half of its heat to its right
neighbor
The Heat Equation: Model 1
• Let Tin denote the temperature at position i at time n.
Particles at the end have a fixed temperature/heat they can
transfer, but they always remain at the same temperature.
n+1
Ti
=
Not transfer heat
0.0Ti−1 + Ti + 0.5Ti+1
n
n
probability 0.75 × 0.5 × 0.25
n
0.5Ti−1
n
n
probability 0.25 × 0.5 × 0.75
+
n
Ti
+ 0.0Ti+1
n
n
n
probability 0.25 × 0.5 × 0.25
n
n
n
probability 0.75 × 0.5 × 0.75
0.5Ti−1 + Ti + 0.5Ti+1
0.0Ti−1 + Ti + 0.0Ti+1
Transfer heat left
n
n
n
probability 0.75 × 0.25 × 0.25
n
n
n
probability 0.25 × 0.25 × 0.75
0.0Ti−1 + 0.5Ti + 0.5Ti+1
0.5Ti−1 + 0.5Ti + 0.0Ti+1
n
0.5Ti
probability 0.25 × 0.25 × 0.25
n
0.5Ti−1
+
+
n
0.5Ti+1
n
n
probability 0.75 × 0.25 × 0.75
0.0Ti−1 + 0.5Ti + 0.5Ti+1
n
n
probability 0.75 × 0.25 × 0.25
n
0.5Ti−1
n
probability 0.25 × 0.25 × 0.75
n
0.0Ti−1 + 0.5Ti + 0.0Ti+1
Transfer heat right
n
+
n
0.5Ti
+ 0.0Ti+1
n
n
n
probability 0.25 × 0.25 × 0.25
n
n
n
probability 0.75 × 0.25 × 0.75
0.5Ti−1 + 0.5Ti + 0.5Ti+1
0.0Ti−1 + 0.5Ti + 0.0Ti+1
The Heat Equation: Model 1
• DEMO
The Heat Equation: Model 2
• Probability is hard!
• It takes a long time
• May only interested in average temperature/heat, not
instantaneous fluctuations
• Take averages
n
Tin+1 = (0.25 × 0.5)Ti−1
+ (0.5 + 0.25 × 0.5 + 0.25 × 0.5)Tin
n
+ (0.25 × 0.5)Ti+1
n
n
Tin+1 = 0.125Ti−1
+ 0.75Tin + 0.125Ti+1
or
Tin+1 = Tin +
n − 2T n + T n
Ti−1
i
i+1
8
• The last line looks like a difference approximation for a
differential equation
The Heat Equation: Model 3
• Let us find a differential equation!
• Make the time increment small
n − 2T n + T n
Ti−1
i
i+1
8
n − 2T n + T n
Ti−1
Tin+1 − Tin
i
i+1
=
δt
8δt
Tin+1 = Tin +
The Heat Equation: Model 3
• Let us find a differential equation!
• Make the space increment small
n
T n − 2Tin + Ti+1
Tin+1 − Tin
= i−1
δt
8δt
Tin+1 − Tin
(δx)2
=
δt
8δt
n −T n
Ti−1
i
δx
−
δx
n
Tin −Ti+1
δx
• Let δx → 0 and δt → 0 such that (δx)2 = δt we get
∂T
1 ∂2T
=
∂t
8 ∂x 2
• Constant in continuum formulation depends on physics
and is usually measured experimentally, or determined
from microscopic (many particle) simulations.
The Heat Equation: Numerical Solution Methods
• Probabalistic
• Finite Difference
• Finite Volume
• Finite Element
• Spectral
Finite Difference Method
• Approximate derivatives by difference quotients
• Simple method to derive and implement
• Convergence rates tend not to be great
• Difficult to use for complicated geometries
• Tends to scale well since communication requirements are
low
Finite Difference Method for Heat Equation
• ut = 8uxx
• Using backward Euler time stepping:
n+1
u n+1 − 2uin+1 + ui+1
uin+1 − uin
= 8 i−1
δt
(δx)2
• Using forward Euler time stepping (strong stability
restrictions):
n
u n − 2uin + ui+1
uin+1 − uin
= 8 i−1
δt
(δx)2
Finite Difference Method for Heat Equation
• Simple method to derive and implement
• Hardest part for implicit schemes is solution of resulting
linear system of equations
• Explicit schemes typically have stability restrictions or can
always be unstable
• Convergence rates tend not to be great – to get an
accurate solution, a large number of grid points are needed
• Difficult to use for complicated geometries
• Tends to scale well since communication requirements are
low
Finite Volume Method for Heat Equation
2
∂u
∂v
•
= ∂∂xu2 or ∂u
∂x = v and ∂t = ∂x
• Consider a cell averaged integral, then use implicit
∂u
∂t
midpoint rule
R xi+1
Rx
ux dx = xi i+1 v dx
Rx
ui+1 − ui = xi i+1 v dx
xi
n+1
n −u n+1 −u n
ui+1
+ui+1
i
i
2
=
n+1
n +v n+1 +v n
vi+1
+vi+1
i
i
δx
2
Rx
ut dx = xi i+1 vx dx
xi
R xi+1
ut dx = vi+1 − vi
xi
R xi+1
n+1
n −u n
ui+1
+uin+1 −ui+1
i
δx
2δt
=
n+1
n −v n+1 −v n
vi+1
+vi+1
i
i
2
• Several ways of approximating the integrals. The one
above is a little unusual, most finite volume schemes use
left sided or right sided approximations.
Finite Volume Method for Heat Equation
• For implicit schemes, hardest part is solving the system of
equations that results
• Explicit schemes parallelize very well, however a large
number of grid points are usually needed to get accurate
results
• Automated construction of simple finite volume schemes is
possible, making them popular in packages
• No convergence theory for high order finite volume
schemes
• Tricky to do complicated geometries accurately
Finite Element Method for Heat Equation
∂2u
• ∂u
∂t = ∂x 2
P
˜ (x, t) = i φi (x)Ti (t), where φi (x) is
• Assume u(x, t) ≈ u
zero for x > xi+1 and x < xi−1 .
• A simple choice is ψ is the triangle hat function, xi+1 − x for
x ∈ [xi , xi+1 ] and x − xi−1 for x ∈ [xi−1 , xi ]
• Now need to find Ti (t), which will be evaluated using finite
differences.
Finite Element Method for Heat Equation
• Consider multiplying the heat equation by a polynomial, φj
that is only non zero over a few grid points, then integrate
by parts,
R xi+1 ∂ 2 u˜
˜
∂u
xi−1 ∂t φi dx = xi−1 ∂x 2 φi dx
R xi+1 ∂ u˜
R xi+1 ∂ u˜ ∂φi
xi−1 ∂t φi dx = − xi−1 ∂x ∂x dx
R xi+1
Finite Element Method for Heat Equation
• Then use implicit midpoint rule
R xi+1
R xi+1 ∂ u˜ ∂φi
= − xi−1
∂x ∂x dx
R xi+1 1 ∂ u˜n ∂ u˜n+1 ∂φi
˜ n+1 −u
˜n
u
φ
dx
=
−
i
δt
∂x + ∂x
∂x dx
xi−1 2
˜
∂u
xi−1 ∂t φi dx
R xi+1
xi−1
1
2
n+1 ˜ n
˜i+1
u
−ui+1
φi |xi+1
δt
n+1 ˜ n
˜i−1
u
−ui−1
φi |xi−1
δt
2δx
∂φi ˜ n+1 ˜n ∂u
1 1 ∂u
+ ∂x 2δx
= −2 2
∂x x
∂x x
xi+1
i+1
i+1
n+1 ∂φi ˜n u
− 21 12 ∂ u˜∂x + ∂∂x
2δx
∂x xi−1
+
xi−1
xi−1
• Since φ are known before hand, can re-write this as matrix
vector products, so need to solve a linear system at each
time step.
Finite Element Method for Heat Equation
• Several other ways of approximating the integrals, can
extend to multiple dimensions.
• Weak formulation allows for solution of equations where
second derivative is not naturally defined
• Large mathematical community developing convergence
theory for these methods
• Well suited to complicated geometries
• Rather difficult to implement compared to other schemes
because of integrals that need to be computed
• Used in many codes, but typically codes are hand written
to obtain high efficiency
• For implicit time discretizations, solving the linear system
of equations that results can be most time consuming part
Fourier Spectral Method for Heat Equation
Fourier Series: Separation of Variables 1
dy
dt
dy
y
Z
dy
y
ln y + a
e
ln y+a
= y
= dt
Z
=
dt
= t +b
= et+b
eln y ea = et eb
eb t
e
y =
ea
y (t) = cet
Fourier Series: Separation of Variables 2
•
∂2u
∂u
=
∂t
∂x 2
• Suppose u = X (x)T (t)
•
dT
dt (t)
T (t)
=
d2 X
(x)
dx 2
X (x)
= −C,
• Solving each of these separately and then using linearity
we get a general solution
•
∞
X
n=0
αn exp(−Cn t) sin(
p
p
Cn x) + βn exp(−Cn t) cos( Cn x)
Fourier Series: Separation of Variables 3
• How do we find a particular solution?
• Suppose u(x, t = 0) = f (x)
• Suppose u(0, t) = u(2π, t) and ux (0, t) = ux (2π, t) then
recall
•
Z
2π
π m=n
,
0 m=
6 n
0
Z 2π
π m=n
,
cos(nx) cos(mx) =
0 m 6= n
0
Z 2π
cos(nx) sin(mx) = 0.
sin(nx) sin(mx) =
0
Fourier Series: Separation of Variables 4
• So if
f (x) =
∞
X
αn sin(nx) + βn cos(nx).
n=0
• then
R 2π
αn =
f (x) sin(nx)dx
R 2π
2
0 sin (nx)dx
0
R 2π
βn =
f (x) cos(nx)dx
.
R 2π
2
0 cos (nx)dx
0
• and
u(x, t) =
∞
X
n=0
exp(−n2 t) [αn sin(nx) + βn cos(nx)]
Fourier Series: Separation of Variables 5
• The Fast Fourier Transform allows one to find good
approximations to αn and βn when the solution is found at
a finite number of evenly spaced grid points
• By rescaling, can consider intervals other than [0, 2π)
• Fourier transform also works on infinite intervals, but
require function to decay to the same constant value at ±∞
Complex Fourier Series
• By using Euler’s formula, one can get a simpler expression
for a Fourier series where sine and cosine are combined
P
• u= ∞
n=−∞ γn exp(inx) x ∈ [0, 2π)
Source: http://en.wikipedia.org/wiki/File:Euler%27s_formula.svg
A Computational Algorithm for Computing An
Approximate Fourier Transform 1
• Analytic method of computing Fourier transform can be
tedious
• Can use quadrature to numerically evaluate Fourier
transforms – O(n2 ) operations
• Gauss and then Cooley and Tukey found O(n log n)
algorithm
• Key observation is to use factorization and recursion
• Modern computers use variants of this idea that are more
suitable for computer hardware where moving data is more
expensive than floating point operations
A Computational Algorithm for Computing An
Approximate Fourier Transform 2
Example pseudo code to compute a radix 2 out of place DFT where x has length that is a power of 2
1: procedure X0,...,N−1 (ditfft2(x,N,s))
2: DFT of (x0 , xs , x2 s, ..., x(N−1)s )
3: if N=1 then
4:
trivial size-1 DFT of base case
5:
X0 ← X0
6: else
7:
DFT of (x0 , x2s , x4s , ...)
8:
X0,...,N/2−1 ← ditfft2(x, N/2, 2s)
9:
DFT of (xs , xs+2s , xs+4s , ...)
10:
XN/2,...,N−1 ← ditfft2(x + s, N/2, 2s)
11:
Combine DFTs of two halves into full DFT
12:
for k = 0 → N/2 − 1 do
13:
t ← Xk
14:
Xk ← t + exp(−2πik /N)Xk +N/2
15:
Xk +N/2 ← t − exp(−2πik/N)Xk +N/2
16:
end for
17: end if
18: end procedure
Sources: http://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm,
http://cnx.org/content/m16336/latest/
The Heat Equation: Finding Derivatives and
Timestepping
• Let
u(x) =
X
ˆk exp(ikx)
u
k
• then
dν u X
ˆk exp(ikx).
=
(ik)ν u
dx ν
• Consider ut = uxx , which is approximated by
ˆk
∂u
ˆk
= (ik)2 u
∂t
ˆkn+1 − u
ˆkn
u
ˆkn+1
= (ik)2 u
δt
ˆkn
ˆkn+1 [1 − δt(ik)2 ] = u
u
ˆkn
u
ˆkn+1 =
.
u
[1 − δt(ik )2 ]
The Heat Equation: Finding Derivatives and
Timestepping
• Python demonstration
The Allen Cahn Equation: Implicit-Explicit Method
• Consider ut = uxx + u − u 3 , which is approximated by
backward Euler for the linear terms and forward Euler for
the nonlinear terms
ˆ
∂u
∂t
ˆ n+1 − u
ˆn
u
δt
ˆ+u
ˆ − uˆ3
= (ik)2 u
ˆ + uˆn − (uˆn )3
= (ik)2 u n+1
The Allen Cahn Equation: Implicit-Explicit Method
• Python demonstration
Key Concepts
• Differential equations used to model physical phenomena -
approximate model
• Microscopic model can suggest ways to find approximate
solutions on computers
• Several methods introduced, probabilistic, finite difference,
finite element, spectral (Fourier transform)
• Need to check simulations for discretization errors, then
before using for prediction, understand modeling errors
• Time stepping, Forward Euler, Backward Euler, Implicit
Midpoint Rule, Implicit-Explicit
References
• Barth, T. and Ohlberger, M. “Finite Volume Methods: Foundation and Analysis”
https://archive.org/details/nasa_techdoc_20030020790
• Boyd, J.P. “Chebyshev and Fourier Spectral Methods” 2nd. edition http:
//www-personal.umich.edu/˜jpboyd/BOOK_Spectral2000.html
• Chen G., Cloutier B., Li N., Muite B.K., Rigge P. and Balakrishnan S., Souza A.,
•
•
•
•
•
•
•
•
West J. “Parallel Spectral Numerical Methods” http://shodor.org/
petascale/materials/UPModules/Parallel_Spectral_Methods/
Corral M. Vector Calculus http://www.mecmath.net/
Cooley J.W. and Tukey J.W. “An algorithm for the machine calculation of complex
Fourier series” Math. Comput. 19, 297–301, (1965)
¨ T. and Herbin, R. “Finite Volume Methods”
Eymard, R. Gallouet,
http://www.cmi.univ-mrs.fr/˜herbin/BOOK/bookevol.pdf
Greenbaum A. and Chartier T.P. Numerical Methods Princeton University Press
(2012)
Heideman M.T., Johnson D.H. and Burrus C.S. “Gauss and the History of the
Fast Fourier Transform” IEEE ASSP Magazine 1, 14–21, (1984)
Ketcheson, D. “Hyperpython” https://github.com/ketch/HyperPython/
Lawler, G.F. “Random Walk and the Heat Equation” American Mathematical
Society (2010)
Macqueron, C. “Computational Fluid Dynamics Modeling of a wood-burning
stove-heated sauna using NIST’s Fire Dynamics Simulator”
http://arxiv.org/abs/1404.6774
References
• McGrattan, K. et al. Fire Dynamics Simulator
http://www.nist.gov/el/fire_research/fds_smokeview.cfm
• Petersen W.P. and Arbenz P. Introduction to Parallel Computing Oxford University
Press (2004)
• Quarteroni, A. Sacco, R. and Saleri, F. “Numerical Mathematics”
http://link.springer.com/book/10.1007%2Fb98885
• Sauer T. Numerical Analysis Pearson (2012)
• Suli,
¨ E. “Finite element methods for partial differential equations”
people.maths.ox.ac.uk/suli/fem.pdf
• Trefethen, L.N. “Finite Difference and Spectral Methods for Ordinary and Partial
Differential Equations”
http://people.maths.ox.ac.uk/trefethen/pdetext.html
• Vainikko E. “Scientific Computing Lectures”
https://courses.cs.ut.ee/2013/scicomp/fall/Main/Lectures
¨
• Vainikko E. Fortran 95 Ja MPI Tartu Ulikool
Kirjastus (2004)
http://kodu.ut.ee/˜eero/PC/F95jaMPI.pdf
Acknowledgements
• Oleg Batra˘sev
• Leonid Dorogin
• Michael Quell
• Eero Vainikko
• The Blue Waters Undergraduate Petascale Education Program
administered by the Shodor foundation
• The Division of Literature, Sciences and Arts at the University of
Michigan
Extra Slides
Fun with Fourier Series
• A Fourier series can represent a regular k-gon in the
complex plane
P
−2 exp [i(1 + kn)t]
• f (t) = +∞
n=−∞ (1 + kn)
t ∈ (−π, π)
Robert, A “Fourier Series of Polygons” Amer. Math. Monthly 101(5) pp. 420-428 (1994)
Schonberg, IJ “The Finite Fourier Series And Elementary Geometry” Amer. Math. Monthly 57(6) pp. 390-404 (1950)
Monte Carlo Method: A Probabilistic Way to Calculate
Integrals
• Recall
¯f =
1
b−a
Z
b
f (x) dx
a
• Hence given ¯
f
Z
b
f (x)dx = (b − a)¯f
a
• Doing the same in 2 dimensions and estimating the error
using the standard deviation
s
ZZ
f (x, y ) dA ≈ A(R)¯f ± A(R)
f 2 − (¯f )2
,
N
R
• Approximate ¯
f by random sampling
¯f ≈
PN
i=1 f (xi , yi )
N
PN
and
f2
≈
2
i=1 (f (xi , yi ))
N
Monte Carlo Method: Python Program
”””
A program t o approximate an i n t e g r a l u s i n g a Monte C a r l o method
T h i s c o u l d be made f a s t e r by u s i n g v e c t o r i z a t i o n , however i t i s
k e p t as s i m p l e as p o s s i b l e f o r c l a r i t y and ease o f t r a n s l a t i o n i n t o
o t h e r languages
”””
import math
import numpy
import t i m e
numpoints =65536 # number o f random sample p o i n t s
I 2 d =0.0
# i n i t i a l i z e value
I2dsquare =0.0
# i n i t i a l i z e to allow f o r c a l c u l a t i o n of variance
f o r n i n xrange ( numpoints ) :
x=numpy . random . u n i f o r m ( )
y =4.0∗numpy . random . u n i f o r m ( )
I 2 d = I 2 d +x∗x +2.0∗ y∗y
I2dsquare = I2dsquare +( x∗x +2.0∗ y∗y)∗∗2
# we s c a l e t h e i n t e g r a l by t h e t o t a l area and d i v i d e by t h e number o f
# p o i n t s used
I 2 d = I 2 d / numpoints
I2dsquare = I2dsquare / numpoints
E s t i m E r r o r =4.0∗numpy . s q r t ( ( I2dsquare − I 2 d ∗∗2)/ numpoints ) # e s t i m a t e d e r r o r
I 2 d = I 2 d ∗4.0
p r i n t ” Value : %f ” %I 2 d
p r i n t ” E r r o r e s t i m a t e : %f ” %E s t i m E r r o r
Listing 1: A Python program which demonstrates how to use the
Monte Carlo method to calculate the volume below z = x 2 + 2y 2 ,
with (x, y ) ∈ (0, 1) × (0, 4).
Sample Results of Monte Carlo Program
N
16
256
4096
65536
∞
Value
37.09
44.49
44.50
43.80
44
Error Estimate
+/- 10.89
+/- 2.40
+/- 0.60
+/- 0.14
0.0
See http://en.wikibooks.org/wiki/Parallel_
Spectral_Numerical_Methods/Introduction_to_
Parallel_Programming for more information