Re: Speeding Up Realtime Visualization
- To: mathgroup at smc.vnet.net
- Subject: [mg84394] Re: Speeding Up Realtime Visualization
- From: Ray Koopman <koopman at sfu.ca>
- Date: Fri, 21 Dec 2007 03:20:52 -0500 (EST)
- References: <fkcugf$5vk$1@smc.vnet.net> <fkeo4m$2gt$1@smc.vnet.net>
On Dec 20, 1:51 pm, Ray Koopman <koop... at sfu.ca> wrote: > 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 ? This may run a little faster. Again, it's untested. While[True, uvv = dt*(u*v^2)[[2;;-2, 2;;-2}]]; {u[[2;;-2, 2;;-2]], v[[2;;-2, 2;;-2]]} = {-uvv + dtf + e*ListConvolve[{{0,1,0},{1,1-4e-dtf,1},{0,1,0}},u], uvv - dtfk + g*ListConvolve[{{0,1,0},{1,1-4g ,1},{0,1,0}},v]}; ugrid = u; vgrid = v ]