Services & ResourcesWolfram ForumsMathGroup Archive

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

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

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

On Sun, 23 Jan 2011, Oliver Ruebenkoenig wrote:

> 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)
>
>
>
>
>

James,

here is a further small speedup:

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 += dt*(du*lap + aIn - u*(k - v*u));
      lap =
       RotateLeft[v, {1, 0}] + RotateLeft[v, {0, 1}] +
        RotateRight[v, {1, 0}] + RotateRight[v, {0, 1}] - 4*v;
      v += dt*(dv*lap + u*(bIn - v*u));
      , {iterationsIn}]; {u, v}],
    CompilationTarget -> "C", RuntimeOptions -> "Speed"
    ];

this changes u = u + ... to u+=... and sets RuntimeOptions->"Speed".


I also tried to parallelize this, but the approach I choose had little 
success:

cflap = Compile[{{uIn, _Real, 1}}, RotateLeft[uIn] + RotateRight[uIn],
    RuntimeAttributes -> {Listable}, CompilationTarget -> "C",
   RuntimeOptions -> "Speed",Parallelization -> True]

re = cflap[u] - 4*u + Transpose[cflap[Transpose[u]]];

I guess here the overhead is too big.

Oliver

• Follow-Ups:
  • Re: Vector problem
    • From: Mark Perrin <m.perrin@me.com>
  • Vector problem
    • From: Mark Perrin <m.perrin@me.com>