CSI 801

Lecture #3

Simulations

September 12, 1996

John Wallin

All that we see or seem
Is but a dream within a dream.

Edgar Allan Poe, A Dream Within a Dream



Pointers and Variables

variables

name associated with a value and a data type

pointers

a name associated with a memory location and a data type


Derived Data Types

data types which are created from several basic or derived data types

allows you to create new data types by encapsulating data

can be made into arrays

can be very useful for science work


Data Structures

dynamically - not predefined and can be updated with new data

organized - created to optimize some set of tasks, usually searching

data sets - groups of basic or derived data types

A way to organize your data


Typical Data Structure Operations

from Cormen et al.

Search(S, k)

Insert(S, x)

Delete(S, x)

Minimum(S)

Maximum(S)

Successor(S,x)

Predecessor(S,x)

where

S = a set of data, k = a key, and x = a data element


Stacks

basic data structures

last in- first out

(LIFO)

input data is pushed on the stack

2 23 56 3

push 5

5 2 23 56 3

output data is popped from the stack

pull x

2 23 56 3

and x=5

similar to HP calculators


Queues

basic data structures

first in- first out

(FIFO)

input data is enqueued

2 23 56 3

enqueue 5

5 2 23 56 3

output data is dequeued

dequeue x

5 2 23 56

and x=3

similar to supermarkets


Link Lists

basic data structures

useful for tracking binned data

---------

1) each bin points to the first element in the bin

2) each element points to the next element in the same bin

3) each element may point to the previous element in the same bin

4) the final element in a bin points to nothing


Binary Trees

basic data structures

1) data is divided into nodes

2) each node has left and right children (pointers)

3) each node has a key used for comparison

4) each node has a parent (pointer)


'The information superhighway will save the environment and revolutionize education in the 21st century."

-An Al Gore -ism


Order of Calculations

How many operations will it take to complete this algorithm?

Binary Tree Searches

log N operations

Finding the minimum of a data set

N operations

Fast Fourier Transforms

N log N

Finding the minimum difference between data values

N^2 operations

Simple inversion of an N x N matrix

N^3 operations


Algorithms

3 key ideas

iteration

recursion

bisection



Bubble Sort Pseudocode

Bubble(A, r)
for i=1, r-1
x = A(i)
for j=i+1, r
if x < A(j) then
exchange A(i) and A(j)
x = A(i)


Bubble Sort in Action

34 35 43 57 12 23 45 23

x = 34

12 35 43 57 34 23 45 23

x=12

12 35 43 57 34 23 45 23

x = 35

12 35 43 57 34 23 45 23

x = 34

12 34 43 57 35 23 45 23

x = 23

12 23 43 57 35 34 45 23

etc...


Quick Sorts

a good algorithm

two parts

a partitioning routine

a quicksort routine

Typical performance

N log N

Worst Case

N^2

example taken from

Cormen, Leiserson, and Rivest

Introduction to Algorithms


Quicksort Pseudocode

QuickSort(A, p, r)
if p < r
 then q = Partition(A,p,r)
Quicksort(A,p,q)
Quicksort(A,q+1,r)

Notes:

initially call with

A = unsorted array

p = 1

r = size of (A)


Partition Pseudocode

Partition(A, p, r)
x = A[p]
i = p - 1
j = r + 1
while TRUE
  do repeat j = j -1
    until A[j] <= x
     repeat i = i + 1
    until A[i] >= x
  if i < j
   then exchange A[i] with A[j]
   else return j
From Cormen et al.

Divides and partitions array into above and below x

Returns the index where this partition occurs


Quick Sort in Action

34 35 43 57 12 23 45 23

x = 34

23 35 43 57 12 23 45 34

23 23 43 57 12 35 45 34

23 23 12 57 43 35 45 34

j=3 is returned

thus

23 23 12 < 34

57 43 35 45 34 >= 34

repeat with smaller arrays


Simulations and Models

Simulation - an attempt to imitate the behavior of a system

Model - the set of dynamical equations used to create a simulation

or visa versa


What are the common characteristics of computational simulations?

Not scientific domain

Not computer algorithm or data structure

Not numerical or mathematical methods


The Gravitational N-body problem

An Example Simulation

The basic scientific question-

How do clusters of stars dynamically evolve?

or

The structure of star clusters does not significantly change with time.


The Gravitational N-body problem

What physical phenomena governs the evolution of this system?

Gravity

Newton's laws

Can this problem be approached analytically?

No. The N-body problem generally cannot be solved?


Mathematical Model

part 1

What is the geometry of the system?

Three dimensional- no symmetry axis


Mathematical Model

part 2


Mathematical Model

part 3

How does this apply to the system?

The forces on a given particle is created by every other particle.

If there are N particles, there are N-1 forces acting on it.


Discrete Algebratic Approximation

We will approximate time as a set of discrete steps.

Vi+1 = Vi + ai+1 dt

Xi+1 = Xi + Vi+1 dt

This is the finite difference system we will use to update positions and velocities.

This is really a numerical solution to the ODE's in the problem.


Simulation Algorithm

What steps will the computer use to solve this problem?

1. Load the particles from a file

2. Calculate the force on each particle

3. Move each particle to a new location

4. Output the result into a data file

5. Repeat 2-4 until the simulation is done.


1. Load the particles from a file

What format is the data going to be in?

How many particles will be allowed in the simulation?

What program generates the initial positions of the particles?

How important are the initial conditions?

Do variables need to be initialized when the program starts?


2. Calculate the force on each particle

for i= 1, n

for j = i+1, n

calculate r(i,j)

ftmp = -G * m(i) * m(j) / r(i,j)*r(i,j)

calculate X(i,j), Y(i,j), Z(i,j)

fX(i) = ftmp * X(i,j) / r(i,j) + fX(i)

fY(i) = ftmp * Y(i,j) / r(i,j) + fY(i)

fZ(i) = ftmp * Z(i,j) / r(i,j) + fZ(i)

fX(j) = -ftmp * X(i,j) / r(i,j) + fX(j)

fY(j) = -ftmp * Y(i,j) / r(i,j) + fY(j)

fZ(j) = -ftmp * Z(i,j) / r(i,j) + fZ(j)


3. Move each particle to a new location

every particle must be moved.

for I = 1, n

VXi+1 = VXi + Fi+1 /mi+1 dt

VYi+1 = VYi + Fi+1 /mi+1 dt

VZi+1 = VZi + Fi+1 /mi+1 dt

Xi+1 = Xi + VXi+1 dt

Yi+1 = Yi + VYi+1 dt

Zi+1 = Zi + VZi+1 dt


4. Output the result into a data file

How often do results need to be saved?

What information is needed in the output file?

How much storage space will be needed for the files?


5. Repeat 2-4 until the simulation is done.

Steps 2-4 must be in a loop.

When is the program done?


Creating the Code

technical plan - part 1

1. Code will be written in ANSI C

2. First targeted platform will be a Sparc 10

3. Simple programming standards will be maintained

4. The code will have internal documentation.

5. The user interface and internal diagnostics will be specified

6. A plan for validation will be created


Code Validation

technical plan - part 2

How do I know this works?

syntax errors

does the darn thing compile?

Fix minor typographical errors

logic errors

are there problems with the coding?

algorithm errors

does the code solve the mathematical model?

model errors

is the mathematical model correct?


Validation and Testing

symmetry

does geometry affect the results

convergence

do the algorithms actually converge

conservation

are conserved quantities conserved?

consistency

is it consistent with previously published results or analytic solutions


Timing and Memory Tests

What is the biggest and bestest simulation you can do? Is this sufficient?

How does the time scale as the problem increases in complexity?

What is the slowest part of the code? Can it be improved?

Is there a memory limitation in the simulation?


Simulation of Phenomenia

Set initial conditions

Run code

Check conservation results

Visualize results


Interpret and Analyze Results

What does it all mean?

Is the hypothesis true or false?

What predictions can you make based on the simulation?

Define and write your conclusions, both technical and scientific!

Remember: We are not just computer geeks, we are also scientists!


Comparing Simulations and Data

Quantitative Comparision

numerical comparison of data and model

how do the values compare?

Qualitative Comparision

less rigourous comparison between data and models

how does the morphology and structure compare?


Models of Stellar Interiors

What is the basic internal structure of a star?


Stellar Structure

What physical phenomena governs the evolution of this system?

Hydrostatic Equilibrium

Mass Continuity

Radiative Transfer

Energy Generation

Equation of State


Mathematical Model

part 1

What is the geometry of the system?

One dimensional- radially symmetric


Mathematical Model

part 2

What equations govern stellar interiors?

A set of five coupled partial differential equations

gravity = pressure

mass is always positive

temperature as a function of radius is related to convection and radiation through the opacity

energy is generated only under certain conditions

opacity is determined by composition, pressure and temperature


Discrete Algebratic Approximation

We will approximate radius as a set of discrete steps.

These may be equally spaced or unequally spaced.

Ri+1 = Ri + dRi

dRi != dRi

PDE's will be solved using first order finite difference solution


Simulation Algorithm

What steps will the computer use to solve this problem?

1. Load initial model

2. Calculate pressure, gravity, temperature, and energy generation

3. Alter star's structure based on gravity and pressure

4. If the structure has signficantly changed, repeat 2 and 3

5. Output final model


1. Load initial model

Guess for the internal mass distribution

Guess the internal temperature profile

Calculate the pressures

Set up reasonable initial conditions and boundary conditions!


2. Calculate pressure, gravity, temperature, and energy generation

Use the discrete algebraic approximation and the mathematical model to determine a new set of physical parameters.

Step through the entire star in concentric shells.


3. Alter star's structure based on gravity and pressure

Use an iterative approach to solving the structure.

4. If the structure has signficantly changed, repeat 2 and 3

Last value creates a new value

This assumes the system converges.

Convergence must be tested.


5. Output final model

Generally, only the final model needs to be saved.

Other iterations are not stable results.

Only one set of data per run is produced.


Technical Plan

Code Validation

symmetry

convergence

conservation

consistency

Timing and memory tests

Simulation

Interpretation and visualization

Analysis and Conclusions


Some final summary points

Computer simulations are NOT just writing some code.

Don't forget the science.

Know thine code.

Remember: We are not just computer geeks, we are also scientists!


Unix This Week

Learn about

diff grep wc cmp cut tr

date du df echo date

ps kill

stty


Projects:

Literature search and definition of your scientific question is due next Thursday.

A few words about plagurism.

Don't do it.

Cite, quote, and give credit to those who deserve it.


Programming Assignment:

random number generators

bubble sort

quick sort

Timing tests:

you do not have to use the exact number of points I specified

Use the 'limit' command to extend your CPU time on the science cluster.

_______________________________

Return to contents page