Simulation of Anomalous Diffusion in Parallel

1. Simulation of Anomalous Diffusion in Parallel Karin Leiderman Final Project – December 3, 2003 Math 471 Dr. Tim Warburton

2. Motivation 1: FRAP Monitoring the fluorescence before photobleaching. Monitoring the fluorescence after photobleaching Monitoring the recovery of fluorescence after photobleaching This is a measurement of the "diffusional mobility" which is usually called lateral mobility since this experiment is almost always done is a planar lipid bilayer.

3. Motivation 2: Single Particle Tracking http://www.jcb.org/cgi/content/full/jcb.200202050/DC1/6

4. Decide on probability distribution function: • The probability of a large jump is low and the probability of a small jump is high • Create a function that returns the jump size when given a random number between 0 and 1.

5. Pseudo Code: • Create a class ‘Particle’ • Holds the x and y coordinates, status and type • Coordinates: x in [0,2n], y in [0,n], n = sqrt(Nprocs) • Status = 0 if the particle will stay on proc and 1 if the particle will leave • Type = 0 if the particle originated outside of the photobleaching region, and 1 if it originated inside the photobleaching region. • Label = particle number • Generate all the particles – equal number on each processor • For each time step { • Run boundaryCheck function: • How many particles are within the maximum jump size of the boundary, or in other words…What is the maximum number of particles that can leave at the next time step? • Generate the new locations • Check how many particles are in the photobleaching region • Run escapeCheck function: • How many particles will actually leave the processor? • Now we can calculate how many particles will come and go from each processor

6. Pseudo Code continued… • For each processor (except self): • As in the sparse matrix-vector multiplication scheme, need to tell each processor how many particles it will send • MPI_Isend this number • MPI_Wait for the send to complete • Now, send the particles…but this is more difficult than it seems • Each vector of particles is of type Particle, so in order to send, we must create our own MPI data type: MPI_Datatype Particles; MPI_Datatype Type[5]={MPI_DOUBLE,MPI_DOUBLE,MPI_INT,MPI_INT,MPI_INT}; int blocklen[5]={1,1,1,1,1}; MPI_Aint disp[5]; int base; MPI_Address(&escapes[0][0].position_x,disp); MPI_Address(&escapes[0][0].position_y,disp+1); MPI_Address(&escapes[0][0].type,disp+2); MPI_Address(&escapes[0][0].label,disp+3); MPI_Address(&escapes[0][0].status,disp+4); base = disp[0]; for (i=0;i<5;i++){ disp[i] -= base; } MPI_Type_struct(5,blocklen,disp,Type,&Particles); MPI_Type_commit(&Particles); class Particle{ public: double position_x; double position_y; int type; int label; int status; };

7. Pseudo Code continued… • For each processor { • MPI_IRecv the number of particles to arrive • Now figure out the new size of the vector of particles • New size = old number of particles + incoming – outgoing • Reallocate space for the new vector, bigger or smaller: Status 0 Status 0 Status 1 Status 0 Status 0 Status 0 for(r=0;r<num_particles;r++){ if(particle[r].status==0){ temp_particle=particle[r]; particle[new_num]=temp_particle; new_num++; } } Status 1 1. Set to a temp 2. Re-fill vector 3. Then realloc Status 0

8. Pseudo Code continued… • For each processor { • MPI_Recv the new particles • Meanwhile… • Each processor is outputting the coordinate information to separate files • Processor 0 holds the photobleaching region and outputs the number of particles inside at each time step Area of photobleaching region / total Area * total number of particles could give a good estimate of the number of particles in the region at any given time This number is ~ 980 particles we can see that this is not the case I will play around with the PDF to see what kind of numbers I can get…

9. Profile: 4 processors, 40,000 particles, 400 time steps Time steps

10. Profile: 16 processors, 40,000 particles, 40 time steps Smaller load on these processors, can we see this from upshot??? Time steps * Number of particles on each processor 96704 proc_0.txt 97101 proc_1.txt 102029 proc_10.txt 101471 proc_11.txt 98779 proc_12.txt 100833 proc_13.txt 101205 proc_14.txt 100352 proc_15.txt 100814 proc_2.txt 99539 proc_3.txt 98584 proc_4.txt 100843 proc_5.txt 101303 proc_6.txt 100740 proc_7.txt 98461 proc_8.txt 101242 proc_9.txt 1600000 total

11. At each time step, each processor has a different number of particles, so the workload is not balanced, but it does not seem to have any significance overall.

12. 400 Particles (gold labeled proteins)40 time steps

13. Just for fun…not an area preserving map

14. Future Directions and Acknowledgements: • I will be working on this code for the next few months as I hope to develop it into s master’s thesis. • As Dr. Warburton pointed out, the need for each processor to have control over each particle and physically send it to another processor is kind of useless for doing the simulation I have just showed you, although this will be very important as I build this code to have clusters and be in 3 dimensions…so all my work is not lost or without purpose!!! • Thanks to Ken Jacobsen’s group at UNC, Chapel Hill, for the movie • Thanks to STMC at UNM Health Sciences Center for the images and data