MathGroup Archive 2011

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

Search the Archive

Re: Do I need MathLink to run finite-difference fast enough for

  • To: mathgroup at smc.vnet.net
  • Subject: [mg115854] Re: Do I need MathLink to run finite-difference fast enough for
  • From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
  • Date: Sun, 23 Jan 2011 05:34:05 -0500 (EST)

On Sat, 22 Jan 2011, James wrote:

> Hello Mike,
>
> This is a Brusselator model of chemical reaction diffussion that is 
> employing the Euler numerical method to solve the PDE system:
>
> u_t=D_u lap(u)+a-(b+1)u+u^2v
>
> v_t=D_v lap(v)+bu-u^2v
>
> For suitable values of the parameters, the characteristic Turing 
> patterns of dots and stripes should emerge.
>
> Using Oliver's Compiler suggestions above under V7, I can run the 
> simulator on a 64X64 grid 10000 times on my machine in about 5.5 
> seconds. 
> However, I feel that's still too slow to create a reasonable interactive 
> Demonstration of the Brusselator where the user will be changing the 
> paramters.  Ideally, I'd like to get it down to about 2 seconds.  I'm 
> not 
> sure how fast it would run in V8.  Here is the complete code using 
> Oliver's suggestions with some additional improvements for stripes that 
> runs in 5.5 seconds on my machine.  Can anyone suggest of a way to get 
> it down to two seconds?
>
> n = 64;
> a = 4.5;
> b = 7.5;
> du = 2;
> dv = 16;
> dt = 0.01;
> totaliter = 10000;
> u = a + 0.3*RandomReal[{-0.5, 0.5}, {n, n}];
> v = b/a + 0.3*RandomReal[{-0.5, 0.5}, {n, n}];
> cf = Compile[{{uIn, _Real, 2}, {vIn, _Real, 2},
>     {aIn, _Real}, {bIn, _Real}, {duIn, _Real},
>     {dvIn, _Real}, {dtIn, _Real}, {iterationsIn,
>      _Integer}}, Block[{u = uIn, v = vIn, lap, dt = dtIn,
>      k = bIn + 1, kern = N[{{0, 1, 0}, {1, -4, 1},
>         {0, 1, 0}}], du = duIn, dv = dvIn},
>     Do[lap = RotateLeft[u, {1, 0}] + RotateLeft[u,
>           {0, 1}] + RotateRight[u, {1, 0}] +
>          RotateRight[u, {0, 1}] - 4*u;
>        u = u + dt*(du*lap + a - u*(k - v*u));
>        lap = RotateLeft[v, {1, 0}] + RotateLeft[v,
>           {0, 1}] + RotateRight[v, {1, 0}] +
>          RotateRight[v, {0, 1}] - 4*v;
>        v = v + dt*(dv*lap + u*(b - v*u)); ,
>       {iterationsIn}]; {u, v}]];
> Timing[c1 = cf[u, v, a, b, du, dv, dt, totaliter]; ]
> ListDensityPlot[c1[[1]]]
>
>

There are two omissions in the above code. a should be aIn and b should 
be bIn.

OK, this fixed code runs on my laptop in V7 in about 9.8 sec.

In V8 this runs in about 6.5 sec. If I then use CompilationTarget -> "C" 
it runs in about 2.9 sec. (actually in a second run it then takes about 
2.6 sec, because the M- and the compiler then have initialized everything 
they need)





  • Prev by Date: Re: Plot works in Mathematca 7 but not in Mathematica 8 [CORRECTION]
  • Next by Date: Re: Converting XML DATEEVENT to Mathematica AbsoluteTime
  • Previous by thread: Re: Do I need MathLink to run finite-difference fast enough for
  • Next by thread: Re: Do I need MathLink to run finite-difference fast enough for