Re: Do I need MathLink to run finite-difference fast enough for
- To: mathgroup at smc.vnet.net
- Subject: [mg115862] Re: Do I need MathLink to run finite-difference fast enough for
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sun, 23 Jan 2011 05:35:37 -0500 (EST)
----- Original Message -----
> From: "James" <icorone at hotmail.com>
> To: mathgroup at smc.vnet.net
> Sent: Saturday, January 22, 2011 2:24:08 AM
> Subject: [mg115846] Re: Do I need MathLink to run finite-difference fast enough for
> 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]]]
Slightly faster, because you sorta left 'a' and 'b' outside:
cf = Compile[{{uIn, _Real, 2}, {vIn, _Real,
2}, {a, _Real}, {b, _Real}, {du, _Real}, {dv, _Real}, {dt, \
_Real}, {iterationsIn, _Integer}}, Module[
{u = uIn, v = vIn, lap, k = b + 1,
kern = N[{{0, 1, 0}, {1, -4, 1}, {0, 1, 0}}]},
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]]]
Two questions:
Do you know why this balks if I raise dt to, say, 0.02? (There may be a sound reason based on PDE stability of this numeric scheme, I'm not sure.) I raise this question because clearly you might be able to get better speed via larger time steps, if the numerics will support that.
(2) Since u is updated and then used in finding the new v, I tried reversing these to see if it would change. It did. Is this expected behavior?
Daniel Lichtblau
Wolfram Research