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}]]