Re: Speeding Up Realtime Visualization
- To: mathgroup at smc.vnet.net
- Subject: [mg84384] Re: Speeding Up Realtime Visualization
- From: Ray Koopman <koopman at sfu.ca>
- Date: Thu, 20 Dec 2007 16:50:33 -0500 (EST)
- References: <fkcugf$5vk$1@smc.vnet.net>
On Dec 19, 9:27 pm, Will <secandum.oc... at gmail.com> wrote: > Since I haven't mentioned this before, I'm trying to make a two species reaction diffusion simulator by approximating the associated PDE using finite difference integration. These 'rough' simulators, e.g.http://texturegarden.com/java/rd/index.html, are typically very fast. I'm going for something with this speed. > > In my code, I've compiled a function which performs the 'next-time-interval' finite difference calculation. Using the Timing function, I found this calculation takes essentially no time. Mathematica usually says 0. sec. However, when I perform this calculation across the entire 100x100 grid, i.e. ~99^2 calcuations, it takes 0.45 sec. > > I originally computed the entire grid using For loops, which required ~0.5-0.6 sec. Figuring that Mathematica doesn't like procedural programming, I tried implementing a Threading technique to apply the finite difference calculation across the grid, which takes 0.45. Of course, this is way to slow for a 'fast' simulator. > > Again, I can speed things up by lowering my grid size, or keeping the grid size and increasing the spatial interval in the integrator, but I'll be losing significant quality. > > I still think the sluggishness has something to do with the way I'm setting up my matrices, which are used to read/write during the finite difference calculation. I figure I'm forgetting something huge > > Ok, here's my code: > Clear["Global`*"]; > N[ > L=100; > ds=1; > length=L/ds; > dt=0.5; > a=0.082;(*u diffusion*) > b=0.041;(*v diffusion*) > f=0.035;(*u feed/rxn rate*) > k=0.0625;(*v death rate*) > dtf=dt f; > dtfk=dt(f+k); > e=a dt/(ds)^2; > g=b dt/(ds)^2; > > uP=Table[1,{i,1,length},{j,1,length}];(*previous u state*) > uC=Table[0,{i,1,length},{j,1,length}];(*current u state*) > vP=Table[0,{i,1,length},{j,1,length}];(*previous v state*) > vC=Table[0,{i,1,length},{j,1,length}];(*current v state*) > c=2;p=1;(*for double buffer*) > u={uP,uC};(*u[[1]] is previous, u[[2]] is current*) > v={vP,vC};(*v[[1]] is previous, v[[2]] is current*) > > (*initial conditions, small square of altered conc values*) > Do[u[[1,i,j]]=0.35,{i,0.4length,0.6length},{j,0.40length,0.6length}]; > Do[v[[1,i,j]]= Sin[0.1i-0.5length]Cos[0.1j-0.5length]^2,{i,1,length},{j,1,length}]; > ]; > > (*finite difference/rxn diffusion calculation*) > Rxn=Compile[{{c,_Integer},{p,_Integer},{i,_Integer},{j,_Integer}}, > u[[c,i,j]]=u[[p,i,j]](1-4*e-dt*(v[[p,i,j]])^2-dtf)+dtf+e*(u[[p,i+1,j]]+u[[p,i-1,j]]+u[[p,i,j+1]]+u[[p,i,j-1]])//N; > > v[[c,i,j]]=v[[p,i,j]](1-4*g+dt*u[[p,i,j]]*v[[p,i,j]])-dtfk+ > g*(v[[p,i+1,j]]+v[[p,i-1,j]]+v[[p,i,j+1]]+v[[p,i,j-1]])//N; > ]; > > (*Rxn Diffusion*) > While[True, > Outer[ > Thread[Rxn[c,p,#1,#2]]&,Range[2,length-1],Range[2,length-1] > ]; > If[p==1,c=1;p=2,c=2;p=1] ; (*c/p inverter*) > ugrid=u[[2]]; > vgrid=v[[2]]; > ] > > Then I Dynamically visualize either ugrid or vgrid using ArrayPlot or Graphics[Raster[]] , which seem take take the same time. > > Any ideas? > > Thanks. Here's some untested code that works with arrays rather than scalars. {p,c} = {1,2}; While[True, u[[c, 2;;-2, 2;;-2]] = u[[p, 2;;-2, 2;;-2]](1 - 4*e - dt* v[[p, 2;;-2, 2;;-2]]^2 - dtf) + dtf + e*ListConvolve[{{0,1,0},{1,0,1},{0,1,0}},u[[p]]]; v[[c, 2;;-2, 2;;-2]] = v[[p, 2;;-2, 2;;-2]](1 - 4*g + dt* u[[p, 2;;-2, 2;;-2]]*v[[p, 2;;-2, 2;;-2]]) - dtfk + g*ListConvolve[{{0,1,0},{1,0,1},{0,1,0}},v[[p]]]; ugrid = u[[c]]; vgrid = v[[c]]; {p,c} = {c,p} ] Or maybe update u & v simultaneously and do away with the double-buffering: While[True, {u[[2;;-2, 2;;-2]], v[[2;;-2, 2;;-2]]} = {u[[2;;-2, 2;;-2]](1 - 4*e - dt*v[[2;;-2, 2;;-2]]^2 - dtf) + dtf + e*ListConvolve[{{0,1,0},{1,0,1},{0,1,0}},u], v[[2;;-2, 2;;-2]](1 - 4*g + dt*u[[2;;-2, 2;;-2]]*v[[2;;-2, 2;;-2]])- dtfk + g*ListConvolve[{{0,1,0},{1,0,1},{0,1,0}},v]}; ugrid = u; vgrid = v ] Is there any particular reason why you don't precompute 1-4*e-dtf and 1-4*g ?