Forside

Abstract

Preface
Beskrive hva prosjektet går ut på. Definere hvor forrige prosjekt sluttet og hvor dette begynner.

Index
Abstract
Preface
Index
1. Introduction
1.1. Previous work
2. The Model
2.1. Basic equations of fluid flow
2.2. Semi-Lagrangian Advection
2.3. The Poisson Equation
2.4. The method of Conjugate Gradients
2.5. GE, LU Factorization
2.6. Bresenham
2.7. Lighting
2.8. Different additional theory goes here
3. Algoritmic and program analysis
3.1. Profiling of 2D prototype
3.1.1. Algorithms
3.1.2. Implementation
3.2. CG benchmarks
3.2.1. Preconditioners for CG
3.3. GE/LU Factorization benchmarks and discussion
3.4. Lighting and raytracing
3.5. Rendering
4. Parallell analysis
4.1. Pseudocode and program analysis
4.1.1. Step 1: Add forces to velocities
4.1.2. Step 2: Solve advection term (Semi-Lagrangian)
4.1.3. Step 3: Solve Poisson Equations
4.1.4. Step 4: Conserve mass by subtracting gradient
4.1.5. Step 5: Advect density and temperature (Semi-Lagrangian)
5. Implementation
6. Results
7. Conclusion
Appendix
References

1.	
Introduction
1.1.	Previous work

2.	
The Model
2.1.	Basic equations of fluid flow
This model assumes that gases can be modeled as inviscid, incompressible, constant density fluids. The effects of viscosity are negligible in gases especially on coarse grids where numerical dissipation dominates physical viscosity and molecular diffusion. When the smoke’s velocity is well below the speed of sound the compressibility effects are negligible as well, and the assumption of incompressibility greatly simplifies the numerical methods. Consequently, the equations that model the smoke’s velocity, denoted by  are given by the incompressible Euler equations. 

Eq 1	

Eq 2	

These two equations state that the velocity should conserve both mass (Eq 1) and momentum (Eq 2). The quantity p is the pressure of the gas and f accounts for external forces. Also the fluid's density is arbitrarily set to one. As in Error! Reference source not found.], Error! Reference source not found.] and Error! Reference source not found.] the equations are solved in two steps. First, an intermediate velocity field u* are computed by solving Eq 2 over a time step t without the pressure term

Eq 3	

After this step, the field u* is forced to be incompressible using a projection method Error! Reference source not found.]. This is equivalent to computing the pressure from the following Poisson equation

Eq 4	
	
with pure Neumann boundary condition, i.e. , at a boundary point with normal n. The intermediate velocity field is then made incompressible by subtracting the gradient of the pressure from it

Eq 5	

Equations for the evolution of both the temperature T and the smoke’s density   are also required. These two scalar quantities are simply moved (advected) along the smoke’s velocity

Eq 6	
Eq 7	

Both the density and the temperature affect the fluid’s velocity. Heavy smoke tends to fall downwards due to gravity while hot gases tend to rise due to buoyancy. A simple model is implemented to account for these effects by defining external forces that are directly proportional to the density and the temperature

Eq 8	

where  points in the upward vertical direction,  is the ambient temperature of the air and  and  are two positive constants with appropriate units such that Eq 8 is physically meaningful. Note that when and , this force is zero. Eq 2, Eq 6 and Eq 7 all contain the advection operator . As in Error! Reference source not found.], this term is solved using a semi-Lagrangian method Error! Reference source not found.], and the Poisson equation (Eq 4) for the pressure is solved using an iterative solver. Section Error! Reference source not found. show how these solvers can also handle bodies immersed in the fluid.

2.2.	Semi-Lagrangian Advection 
2.3.	The Poisson Equation
2.4.	The method of Conjugate Gradients
2.5.	GE, LU Factorization
2.6.	Bresenham
2.7.	Lighting
2.8.	Different additional theory goes here

3.	
Algoritmic and program analysis
3.1.	Profiling of 2D prototype
3.1.1.	Algorithms
Presentere benchmarkresultater fra prototypen. Vise hvilke områder som kan/bør optimaliseres og diskutere hvordan det evt kan gjøres. 
3.1.2.	Implementation
Analysere implementeringsmessige løsninger i prototypen, og hvordan disse påvirker hastigheten. Diskutere forbedringer og nye elementer ved utvidelse til tre dimensjoner. Analysere hvilke data modellen krever, og etter hvordan mønster de vil aksesseres. Videre diskutere forskjellige minnestrukturer med hensyn på hvordan de vil prestere hastighetsmessig (hovedsaklig mhp cache).  
3.2.	CG benchmarks
Fordeler, ulemper, diskusjon. 
3.2.1.	Preconditioners for CG
3.2.1.1.	PETSc Preconditioners
PETSc implements a number of different preconditioners for CG, together with a few direct solvers. Table 1 shows a comparison of the different PETSc preconditioners. Two of the most promising preconditioners (Incomplete Cholesky and Multigrid) aren't implemented for the testprogram's matrix format, and hence not included in the table. See section XX for a description of these two. 
Preconditioner	Time	Iterations	Error	
None	13.25 s	475	0.000944	
Jacobi	13.59 s	475	0.000944	
SOR	38.92 s	1000	34101857.585735	
Block Jacobi	7.40 s	142	0.000690	
Eisenstat	11.49 s	168	0.000747	
Redundant	8.66 s	142	0.000690	
Asm	9.89 s	142	0.000690	
Incomplete LU	8.50 s	142	0.000690	
Table 1 Comparison of PETSc's preconditioners

Direct solvers				
SLES	55.49 s	2	0.000000	
LU	4.75 s	1	0.000000	
Table 2. 

The PETSc solver was tested with a set of equations with 40.000 unknowns, built out of a standard 2D Laplace problem. All benchmarks were run with a suggested solution vector. This vector was exported from the prototype, hence representing an actual situation in the simulator.

The two 

Fordeler, ulemper, diskusjon. Presentere de forskjellige som er testet i PETSc, og resultatene derfra.  
3.3.	GE/LU Factorization benchmarks and discussion
Fordeler, ulemper, diskusjon. Presentere resultater fra testprogram, minnebruk, prosesseringstid etc etc. 
3.4.	Lighting and raytracing
Analysere dataaksessmønster og forventet prosesseringstid ved forskjellige lysmodeller, forskjellige lys, varierende antall lys. Forskjellige linjealgoritmer (f.eks. Bresenham). 
3.5.	Rendering
Metoder og algoritmer for rendering. Analysere tidsforbruk ved sending av datamengder over minnebuss, forskjellige OpenGL-løsninger.
4.	
Parallel analysis
4.1.	Parallel system
TODO: Describe the type of parallel system we're designing for here.... One process, several threads, several cpus, shared memory etc etc.
4.2.	The data
TODO: Describe how the prototype's data were organized. Give a brief description of the model's 3D version. Needed for getting the following sections into their appropriate context. 
4.3.	Pseudocode and program analysis
Top level pseudocode for the entire program is presented in  Figure 1. 
For each frame
1	Add forces to velocities
2	Solve advection term (Semi-Lagrangian)
3	Solve Poisson equations
4	Conserve mass by subtracting gradient
5	Advect density and temperature (Semi-Lagrangian)
6	Calculate lighting
7	Render
Figure 1 Top level pseudocode.
Each of the seven steps are specified in the following sections. 
4.3.1.	Step 1: Add forces to velocities
Step 1 updates the velocity field, according to different forces. This is done simply by multiplying the force of each voxel with the timestep t and adding the result to the voxel's velocity (see Appendix A). The forces include user defined fields and the buoyancy force defined by Eq 8.

For each voxel
1.1	
1.2	
1.3	
Figure 2 Low level pseudocode for Step 1
In other words, a standard vector-vector addition. A simple data-parallel approach would do for this step of the algorithm. Some care should be taken to avoid several threads trying to access the same cacheline at the same time. 
4.3.2.	Step 2: Solve advection term (Semi-Lagrangian)
This step solves for the advection term in Eq 3, using a Semi-Lagrangian Semi-Implicit (SLSI) scheme. This builds a new grid of velocities from the ones already computed, using the trapezoidal method as integration method. The center of each voxel is traced through the velocity field, and the velocities  are linearly interpolated at the arrival points. The values are then transferred back to the cells the trajectories originated from. Simple linear interpolation is easy to implement, gives satisfactory results and is unconditionally stable as it never overshoots data. Higher order interpolation schemes are, however, desirable in some cases for high quality animations, but as this prototype aims for a realtime environment, the processing cost of higher order interpolation is considered not worth the effort. 

For each voxel at position 
2.1	Calculate  
2.2	Set temporary velocity  for each voxel to 

Figure 3 Low level pseudocode for Step 2

Step 2.1 involves linear interpolations of two velocity vectors. The first, , is the current velocity at the voxel to be updated. The other, , is an off-grid velocity vector needed in the trapezoidal method. At first glance, step 2 seems just as suitable for a simple data-parallel approach as step 1. Unfortunately, the off-grid interpolation involves several other voxels in addition to the current voxel. Depending on the length of the velocity vector, the additional voxels need not even be neighbouring voxels but could be located far away from the voxel to be updated. If a data-parallel approach is chosen, this fact may cause race conditions as several threads may try to access a given voxel at the same time.  Restricting the user defined velocities/forces would reduce the possibility of this condition to occur, but such limitations would also reduce the simulator's useability and hence isn't very desirable. In general, these race conditions may be reduced by giving the different threads working areas located sufficiently far away from each other in memory. The distance between the memory locations depends on the relation between total number of threads and number of voxels in the grid. More threads/less voxels will result in smaller working areas, thus locating them closer to each other and higher possibility of race conditions. 
4.3.3.	Step 3: Solve Poisson Equations
In step 4, the velocity field is forced to conserve mass, and this requires the solution of a Poisson equation for the pressure (Eq 4) in Step 3. A matrix of coefficients is constructed using the five-point formula approximation to Laplace's equation, and the system of equations is solved using the Conjugate Gradient method (CG). 

3.1	Build solution vector . For each voxel 
: 
3.2	Solve equations (CG):

Figure 4 Low level pseudocode for step 3
4.3.4.	Step 4: Conserve mass by subtracting gradient
Step 4 is a standard vector-vector subtraction:

For each voxel 
4.1 Update velocity 
4.3.5.	Step 5: Advect density and temperature (Semi-Lagrangian)
For each voxel at position 
5.1	Calculate  
5.2	Set temporary velocity  for each voxel to 
Figure 5 Low level pseudocode for step 5

5.	
Implementation
Beskrive alle relevante aspekter ved implementasjonen som ikke allerede har kommet frem i tidligere seksjoner.
6.	
Results
Visuelle resultater. Ytelsesmessige resultater. 

7.	
Conclusion

Appendix

References