Transcript Numerical Simulations

Swinburne Astronomy Online
Module:
Computers and Astrophysics
Activity:
Astrophysical Simulations
© Swinburne University of Technology
Summary
In this Activity, we will learn about numerical simulations
and their importance in astrophysics today. In particular,
we will look at:
(a) why we use numerical simulations in astrophysics
(b) how to set up a numerical experiment
(c) N-body codes for dynamic simulations
(d) examples of N-body simulations
(e) adding gas
Introduction
Computer experiments - or numerical experiments - are
mathematical models for simulating the behaviour of
complex physical systems. There are various reasons
why one may want to use a computer experiment,
including:
• to help design a complex device before using expensive
technology to actually build it;
• to gain information about physical systems that are not
readily available in laboratory experiments (astronomy is
a perfect example of this!); and
• to help understand extremely complex system
(astronomy is also a good example of this).
Dynamic systems that evolve with time, such as
galaxies, have three coordinates of position (since they
are 3D objects) and also three coordinates of velocity
(since they can move in all 3
directions).
However, we can only observe
two position coordinates and the
line-of-sight velocity.
observer
So to determine the underlying structure of the Universe,
we need to use computer simulations.
There are many situations in nature where we can view a
physical system as a collections of particles that interact
and obey certain laws.
Astronomy offers many excellent example of this: the
planets in the Solar System, the stars in the galaxy, and
galaxies in the Universe obeying Newton’s law of
gravitation.
These particles have various physical attributes, such as
mass, position and velocity, and each particle interacts
with all others under the force of gravity. This is known
as the N-body problem (where N is the number of
particles used in the computer experiment), and we shall
explore the N-body problem in this Activity.
Setting up a Numerical Experiment
The general steps involved in running your computer
experiment are as follows:
Physical phenomena
chose the physical system
that you wish to investigate
Setting up a Numerical Experiment
The general steps involved in running your computer
experiment are as follows:
Physical phenomena
Mathematical model
the physical system can be
approximated by a
mathematical model, which
generally uses some simplifying
assumptions to describe the
workings of the physical system
Setting up a Numerical Experiment
The general steps involved in running your computer
experiment are as follows:
Physical phenomena
Mathematical model
Discretise the model
the mathematical model must
be converted from a continuous
or differential equation into an
algebraic approximation which
the computer can solve. One
must be careful at this step to
introduce as few errors as
possible. Both time and space
must be discretised, which can
produce huge numerical errors
Setting up a Numerical Experiment
The general steps involved in running your computer
experiment are as follows:
Physical phenomena
Mathematical model
Discretise the model
Numerical algorithm
chose the choice of
discretisation is often related
to the algorithm chosen to
solve the discrete system.
You need to be able to solve
the discrete problem rapidly
on your computer, else it’s
all useless!
Setting up a Numerical Experiment
The general steps involved in running your computer
experiment are as follows:
Physical phenomena
Mathematical model
Discretise the model
Numerical algorithm
Computer program
writing the code is what it’s all
about! The code needs to be
well engineered to make use of
the available computing power,
be easy of use and easy to
modify
Setting up a Numerical Experiment
The general steps involved in running your computer
experiment are as follows:
Physical phenomena
Mathematical model
Discretise the model
Numerical algorithm
Computer program
Computer experiment
finally you get to run your
computer experiments, but you
writingtothe
need
know
code
what
is what
you’re
it’s all
about! for,
testing
Thewhat
codeyou’re
needstrying
to beto
well engineered
find,
and how to to
domake
analysis
useon
of
the data
available
that your
computing
experiment
power,
be easy of use
produces.
Pretty
and
pictures
easy toare of
modify vital at this stage!
course
Discretising the Mathematical Model
A major task when setting up a computer experiment is
“discretising” the mathematical model.
Mathematics usually involves continuous functions, whereas
computers can only work with individual bits of information.
For example, a computer does not really have a concept of a
smooth sine wave:
But it can handle a set of discrete points that represent the
continuous sine wave.
The computer can then deal with these discrete bits of
information, adding them, multiply them, storing them etc.
Dynamical Simulations
While there are many systems in astronomy which we
can model using a computer code, in this Activity we will
be concerned with dynamic simulations, which means
the study of systems that evolve with time. We will not
cover static simulations, which study a specific state of a
system.
In most astrophysical systems, the global properties of
the system are determined by the force of gravity. While
there are also other forces involved - like electrical and
chemical forces holding atoms and molecules together it is gravity that is the dominant long range force in the
Universe.
Examples in Astronomy
Dynamical systems in astronomy that are governed by gravity
that we can model using numerical simulations include:
• Solar System studies:
–
–
–
–
the study of the dynamics and stability of the Solar System
the motions and interactions of the asteroids in the asteroid belt
modelling planetary rings
the orbits of comets
• Star systems:
– the formation of stars through the collapse of
giant molecular clouds
– the dynamics of binary systems and small
stellar clusters
– the formation and evolution of globular clusters
Globular cluster M15
• Galaxy studies:
–
–
–
–
the dynamics of spiral galaxies
interactions between galaxies and galaxy mergers
the evolution of galaxies
galactic clusters and
cluster dynamics
• Cosmological studies:
– the formation of structure
in the Universe
– the dark matter problem
– the evolution of the Universe
Merging galaxies NGC2207 and IC2163
These are just some of the systems that have been modelled
over the past 20 years using dynamic computer simulation. So
let’s have a look at how dynamic computer codes work.
The Mathematical Model
There are two main parts to these dynamical codes:
• the force calculation and
• the time evolution.
Both of these components can be described by a
mathematical model. This means that we have a set of
mathematical equations to solve which will tell us about the
future state of our system given a set of initial conditions.
Once we have our mathematical model, we need to
discretise it and use a numerical algorithm to solve the
discrete system on a computer.
Let’s have a look at the two parts of the mathematical
model that describe a dynamic gravitational system.
Gravitational Forces
Newton’s universal law of gravity for the force between
two bodies is given by
Gm m
1
F=
2
r2
where F is the resulting gravitational force,
G is the universal gravitational constant,
m1 and m2 are the masses of the two bodies, and
r is the separation between the bodies.
m1
r
F1
F2
m2
Forces F1 and F2 are equal in size and opposite in direction.
What if we have five masses in our system?
m4
m3
F13
m2
F14
F12
m1
F15
m5
The force exerted on body 1 by the other 4 bodies would
be given by:
F1 =
G m1 m2
r122
+
G m1 m3
r132
+
G m1 m4
r142
+
G m1 m5
r152
We also need to calculate the force on particle 2 due to the
other 4 particles:
and the force on particle 3 due to
all the other particles:
and the force on particle 4 due
to all the other particles:
m4
m3
m2
m1
m5
and the force on particle 5 due to all the other particles:
You can see that things get quite complicated very quickly
and a computer could be of use!
We can write the force equation quite simply as:
N
Fi = j=1

G mi mj
rij2
where the  (“sigma”) symbol means to sum or add over all
other particles.
We also use subscripts i and j to represent the particles. So
for each N particle i we need to sum over all the other N-1
particles, which we assign with a j subscript.
You can see that the mathematical model for gravitational
force is quite easy to discretise for our N particles. There
are a couple of slight complications however...
Firstly, force is actually a vector quantity, not a scalar
quantity, which means it has a magnitude and a direction.
Our current equation for the force does
N
G mi mj rij
F
=
-
not tell us anything about the direction of
i
j=1
| rij3 |
the force, so that needs to be amended.
Secondly, as the particles get closer together, the forces get
larger. In fact, as particle i and j approach each other the
denominator rij of the force equation approaches zero and so
the force become infinite!
Now computers do not like having to calculate infinite
numbers and so we need to soften the gravity. Our equation
becomes:
N
Fi = - j=1

G mi mj rij
| rij3 | + 
This might not look like much of a change (since rij/|rij|3 ~ 1/rij2),
but it is! The funny vertical bar symbols |rij| are called “mod”
or “absolute value” and they take only the positive value, for
example |-3| = 3.
That means that the remaining rij gives the direction of
the force. If we imagine that the position rij is in 2D
space, with an x-direction
y
x component
component and a y-direction
component, then rij will give
y component
rij
the force its direction.
x
The softening parameter  has to be carefully chosen - if
it’s too large it will affect the physics of the system, and if it
is too small then forces will become huge.
From the list of astrophysical examples of N-body system
we gave earlier, you should be able to start seeing how
each object (planet, star or galaxy) can act as a single
“particle” with its own mass and position. And we’ve just
seen how to calculate its force.
In the case of the Solar
System, we would only
need the Sun and the
nine planets - so N=10.
To be more accurate, we could include the satellites, so let’s
say N=80. And what about the asteroids? We might want
to include 100,000 of them, so now N=100,080. You can
see that the number of bodies in our system can rapidly
increase!
If you want to model a small cluster of stars, you might
want to use a few hundred or few thousand stars. For
globular clusters you will need a few hundred thousand
particles.
This means that the force on each star in the globular
cluster depends on the position and mass of the other
hundred thousand stars in the cluster. You can see that
a lot of calculating needs to be done! In fact, there will
be N x (N-1) calculations, or O(N2) (pronounced “order N
squared”).
Sometimes, however, you just cannot use a particle for
every object in your system. You will not be able use a
hundred million particles to represent all the stars in the
Galaxy. So sometimes your “particles” actually represent
a collection of physical bodies.
The Equation of Motion
The evolution of the system is governed by the equations
of motion, which tells us how the system evolves with time.
Acceleration is defined as the rate of change of velocity,
and velocity is the rate of change of position as follows:
dv
a
dt
dr
and v 
dt
where a = acceleration v = velocity
r = position
t = time
We can easily discretise this mathematical model by:
v vf - vi
r rf - ri
a
and v 


t t f - t i
t t f - t i
where
i=initial f=final
The  symbol (“delta”) represents a small but finite change.
Newton’s second law relates force to motion via the
equation F = m x a. We can use Newton’s law of
gravitation to replace F by Fgrav, which leads to:
dv
F  ma  m
dt
dv Fgrav
Gm 2

- 2
dt
m
r
We want to solve for both the position and velocity of the
system at the next timestep.
F v
Ft
 v 
We can do this by solving the a  
m t
m
following equations:
r
v
 r  vt
t
If we now assume that we take a set timestep t between
the old and new states of the system, and let v = vf - vi
(where vf is the final velocity and vi the initial velocity) and do the
same thing for the position, then
Ft
the new velocity and position is
vf  vi +
given by:
m
rf  ri + v i t
So once we have solved for the gravitational force F at
the initial state of the system, we can solve for the
system at the next timestep.
There are actually many fancier ways to discretise the equations of
motion to give more accurate time stepping, but we won’t go into
them just here.
Accuracy and Stability
There are two more things that we need to be careful
about: our choice of t and N.
The timestep t controls the stability of the simulation. If
we use a constant timestep t in the equations of motion,
we will get large errors when two particles get close. So
we need to use a numerical scheme with a variable
timestep which will automatically decrease t if particles
are too close and increase t as particles move apart.
The particle number N gives the resolution of the
simulation. Obviously we will need more that a few
hundred particles to accurately resolve the dynamics of a
galaxy, but of course there are limits as to how large N
can be. Supercomputers can help us out here.
The N-body Algorithm
So we’re now armed with our mathematical model for the
gravitational force and equations of motion; we have a discrete
algorithm for the mathematical model, and we’re ready to write
our computer code to run our computer experiments.
This is what our computer algorithm will look like:
Set initial conditions
Choose N and t, set initial
particle mi, ri, vi, Fi
Solve equations of motion
ai = vi/ ti & vi = ri/ ti
Calculate forces
Fi = j Gmimj/rij2
Update time counter
tnew = told + t
Output data
rnew, vnew, Fnew, tnew
Examples of N-body Simulations
In the next few slides, we will present some examples of
N-body codes that have been used to simulate the
dynamics of spiral galaxies, the formation of structure
within the Universe, galaxy clustering and the formation
of quasars.
The number of particles in the simulation, N, provides the
resolution, and increasing N allows more detail to be
seen in the simulations.
When N is much larger than about 30,000, we need a
supercomputer for its calculational speed, plus rapid
storage of and access to the data.
Galaxy Dynamics
This is a low resolution (N=7000)
simulation of a spiral galaxy. The
simulation began with a smooth
distribution of particles, and the
action of gravity sets up spiral
arms within the galaxy. Colour
represents density.
Simulations like this have
confirmed theories of spiral arm
structures within galaxies and
answered questions about the
missing mass in the Universe.
Dark Matter in the Universe
This is a low resolution
(N=32,768) simulation of
the dark matter distribution
in a box with a side length
of 100 Mpc*, or 300 million
light years. Each particle
represents a “blob” of dark
matter with a mass ten
times that of the Milky Way.
What are most interesting to
note are the filamentary
structures connecting
regions where the particle
density is much higher than
the average density within
the box. This box does not
represent the entire
Universe!
* Mpc = mega parsec = 106 pc
Click here for animation
Particle Number
These images show the end state of a simulation of a cluster of galaxies
comparable to the Virgo cluster. The colour represents density and the radius of
the computational box is 4 Mpc.
The image on the left is a low resolution (N=67,000) CRAY simulation, while the
image on the right is from a high resolution (N=1.3 million!) parallel simulation
run on a “cluster” supercomputer. As you can see, as more particles are used,
more detail of the cluster’s structure can be seen.
Quasar Formation
in the Early Universe
One of the current problems in
cosmology is this: simulations
have shown that while structure
forms very early in the Universe,
galaxies tend to form quite late.
Yet we know that quasars exist at
very high redshifts. So how did
they form?
This sequence of images shows
the formation of a quasar at a
redshift of z=8 (which means very
early in the history of the
Universe) using a high resolution
simulation.
Improving the Algorithms
We’ve only considered N-body algorithms here, but there other types of
algorithms currently used in astrophysics for solving for the gravitational
force in a physical system.
None of the simulations you have just seen actually use the “direct
summation” technique that you have learned about in this Activity, where
each of the N particles calculates its force due to the N-1 other particles.
This requires N2 calculations, and for N=5000, there are 25 million
calculations to do, and for N=50,000, 2.5 billion calculations are needed!
It just takes too long to run a simulation using this direct approach. And
that’s just the gravitational part of the algorithm: we also need to
integrate the equations of motion. If you have 50,000 particles coming
close together you have to decrease the timestep t.
So computational astrophysicists spend a lot of their time writing very
clever codes to calculate the gravity as quickly and efficiently as
possible. This is usually done by approximating the forces from particles
far away and accurately calculating (by direct summation) the forces of
nearby particles.
Adding Gas
We have looked at several examples of astrophysical
N-body systems where we calculate the force on a set of
gravitating particles that act as point masses.
There are many cases in astronomy where we would also like
to run computer experiments of the gas behaviour of systems.
Examples include the collapse of molecular cores to form
protostars and their surrounding disks; the modelling of gas
jets emanating from young protostars; the expansion of shock
waves into the interstellar medium from supernova explosions;
and the study of the gas disks around cataclysmic variables
and black holes.
We might also want to add gas to our simulation of galaxies,
so that we can see how the stars and gas in the arms of
spirals interact.
The study of the motion of a fluid is called hydrodynamics.
Fluids - both liquids and gases - are governed by pressure
and viscous forces.
Purely pressure forces are expansive forces, compared
with the inward force of gravity. Viscosity, which depends
on the properties of the gas, acts as a type of resistance
in a fluid. Purely viscous forces are spreading forces.
The total force on a particle of gas is given by:
Ftotal = Fgravity + Fpressure + Fviscous
This means that we need to include more physics into
our computer code. We need a mathematical model for
both pressure and viscosity, which we can then
discretise.
Each particle in a gas simulation represents a collection
of molecules or atoms, since of course we cannot
represent each individual gas molecule, just like we can’t
represent each individual star in a galaxy!
In gas simulations, each particles represent a sphere of
gas, which we can think of as a having a “radius of
influence”. Any other particle in particle i’s sphere of
influence is included in all force calculations with it.
Those outside the sphere of influence are not.
interaction radius
i
not used to calculate
the force on particle i
not used to calculate
the force on particle i
Massive Protostellar Disks
This is a simulation of a massive
protostellar disk, which forms
around a young stars as it
collapses out of a molecular
cloud.
The disk is coloured by density
and has the same mass as the
central protostar. The gravity
within the disk causes
instabilities which create
transient spiral arm features and
small blobs to break off and
merge. Could these blobs be
forming planets?
Click here for animation
Disks Around Young Binary Stars
This next simulation shows the
dynamics of gas in a
circumbinary disk, which is a
disk that surrounds a young
binary star system.
At the centre there are two stars
in mutual orbit around each
other, and the action of gravity
sweeps out a hole in the disk,
cutting it off from the central
binary. You can see transient
spiral arms occasionally feeding
the young stars with gas.
Click here for animation
Summary
In this last Activity of the Unit we have seen how
computer experiments are a fantastic tool in modern
astronomy, allowing us to understand the structure and
dynamics of many systems in the Universe.
We have learned about how to set up an N-body
computer experiment and seen several examples of their
use in astronomy today.
With the help of supercomputers, we can now increase
the value of N into the millions, allowing the Universe to
be studied in great detail. It is up the computational
astrophysicists to write better algorithms and codes to
make use of this great computing power!
Image Credits
Sunflower galaxy M63 - Satoshi Miyazaki, Suprime-Cam, Subaru Telescope, NOAJ
http://antwrp.gsfc.nasa.gov/apod/ap000627.html
Globular cluster M15 - Hubble Heritage
http://heritage.stsci.edu/public/2000aug3/ngc7078.html
Interacting galaxies NGC2207 and IC2163 - Hubble Heritage
http://heritage.stsci.edu/public/99nov4/ngc2207.html
The planets - Welcome to the Planets, California Institute of Technology
http://pds.jpl.nasa.gov/planets/
Galaxy simulation - Sarah Maddison, Swinburne University of Technology
(volume rendering by Robin Humble, CITA)
Image Credits
Cosmological simulation and animation - by Chris Fluke of Swinburne University,
using the Hydra N-body SPH code of H. M. P. Couchman & P. A. Thomas
http://hydra.mcmaster.ca/hydra/index.html
Clustering simulation - HPCC, University of Washington
http://www-hpcc.astro.washington.edu/TSEGA/picture/picture.html
Quasar formation simulation - HPCC, University of Washington
http://www-hpcc.astro.washington.edu/TSEGA/picture/picture.html
Massive protostellar disk simulation - Sarah Maddison, Swinburne University of
Technology
Binary disk simulation - Sarah Maddison, Swinburne University of Technology
End of Activity and Unit