Re: NDSolve memory management problem
- To: mathgroup at smc.vnet.net
- Subject: [mg102358] Re: NDSolve memory management problem
- From: schochet123 <schochet123 at gmail.com>
- Date: Thu, 6 Aug 2009 06:33:38 -0400 (EDT)
- References: <h41f71$rjo$1@smc.vnet.net>
On Jul 20, 1:00 pm, David Szekely <dr.szek... at gmail.com> wrote:
> I'm currently trying to run a simulation in Mathematica usingNDSolvein se=
veral >dimensions. However, it tends to chew up systemmemoryvery quickly an=
d give the >standard:
> "No morememoryavailable.
> Mathematica kernel has shut down.
> Try quitting other applications and then retry."
> message.
> Here's some code for the heat equation in 4 dimensions over relatively sm=
all >ranges to show what I mean:
> i = 0;
> sol =NDSolve[{D[u[t, x, y, z], t] ==
> D[u[t, x, y, z], x, x] + D[u[t, x, y, z], y, y] +
> D[u[t, x, y, z], z, z],
> u[0, x, y, z] == 0,
> u[t, 0, y, z] == Sin[t],
> u[t, 40, y, z] == 0,
> u[t, x, 0, z] == Sin[t],
> u[t, x, 40, z] == 0,
> u[t, x, y, 0] == Sin[t],
> u[t, x, y, 40] == 0},
> {u}, {t, 0, 100}, {x, 0, 40}, {y, 0, 40}, {z, 0, 40},
> MaxSteps -> Infinity, MaxStepSize -> 1,
>
> EvaluationMonitor :> If[t > i,
> Print[{t, MemoryInUse[]/1024^2 // N}];
> i += 10;]
>
> ]
> MemoryInUse[]/1024^2 // N
>
> The only solutions I need are at t = 0, 5, 10, 15 etc.. Ideally, what i=
'd like to be >able to do is add a function in "EvaluationMonitor" which re=
moves all values of t >other than the ones that I need *as it solves the sy=
stem*. This would significantly >reduce the strain on systemmemory. I would=
greatly appreciate help on this >problem! Adding morememoryis not a good s=
olution since my actual problem kills >a system with 12 Gb ofmemory:(
There are two steps that have to be taken to get around this problem.
The first is to make NDSolve stop whenever a lot of memory is being
used. The way to do that is with the EventLocator method. The second
step is to make Mathematica give back most of the memory it is using,
while keeping the computed value of the solution. That is slightly
tricky, since if we define something like solnow[x_,y_,z_]=uu
[tnow,x,y,z] where uu is a name given to the InterpolatingFunction
returned by NDSolve then solnow uses the same amount of memory as uu
so doing Clear[uu] won't reduce memory usage. It is therefore
necessary to extract the grid points and values from the
InterpolatingFunction at the current time and use Interpolation to
produce an InterpolatingFunction using only the values at the current
time. The modification of your example given below yields the solution
(called solnow[x,y,z]) at the final time without running out of
memory. You will need to modify it if you also want to store the
solution at particular intermediate times.
Although not strictly necessary, I replaced EvaluationMonitor with
StepMonitor to reduce the number of evaluations; for example, a fourth-
order RungaKutta method evaluates the functions appearing in the ODE
four times for each time step.
In a more complicated problem it may also be necessary to use a 64 bit
system to avoid kernel memory limits, since the size of the steps
needed to keep memory usage down might otherwise be too small.
Presumably the OP uses such a system since he has 12GB memory.
Steve
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]
i = 0; maxmemorytouse = 300 1024^2;
timenow = 0; maxtime = 100; solnow = (0 &);
While[timenow < maxtime,
usol = u /.
NDSolve[{D[u[t, x, y, z], t] ==
D[u[t, x, y, z], x, x] + D[u[t, x, y, z], y, y] +
D[u[t, x, y, z], z, z],
u[timenow, x, y, z] == solnow[x, y, z], u[t, 0, y, z] == Sin[=
t],
u[t, 40, y, z] == 0, u[t, x, 0, z] == Sin[t],
u[t, x, 40, z] == 0, u[t, x, y, 0] == Sin[t],
u[t, x, y, 40] == 0}, {u}, {t, timenow, 100}, {x, 0, 40}, {y, 0,
40}, {z, 0, 40},
Method -> {"EventLocator",
"Event" :> MemoryInUse[] > maxmemorytouse},
MaxSteps -> Infinity, MaxStepSize -> 1,
StepMonitor :>
If[t > i,
Print[{"StepMonitor: At time ", t, "Memory in use is",
MemoryInUse[]/1024^2 // N}];
i += 10;]][[1]];
Print[{"At time ",
timenow =
InterpolatingFunctionDomain[First[Flatten[{usol}]]][[1, -1]],
"Memory in use was ", MemoryInUse[]/1024^2 // N}];
xyzvals = Rest /@ Flatten[Last[InterpolatingFunctionGrid[usol]], 2];
uvals = Flatten[Last[InterpolatingFunctionValuesOnGrid[usol]]];
Clear[usol]; solnow = Interpolation[Transpose[{xyzvals, uvals}]];
Print[{"Memory in use reduced to ", MemoryInUse[]/1024^2 // N}]]