Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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 ]


  • Prev by Date: Re: FindRoot / Jacobian
  • Next by Date: Re: using Dsolve for equation: y'' + |y|^2y=ay
  • Previous by thread: Re: Speeding Up Realtime Visualization
  • Next by thread: packages not working under Linux