How to use NDSolve with limited RAM?

Question

How can I use NDSolve when running out of RAM?

The code I am implementing:

sol = NDSolve[ pdes, {A, B, C},
    {x, -L, L}, {t, Tf, Tf} ]

in order to obtain the end solution of NDSolve to be exported to the next run (due to limited RAM).

In previous attempts the solution of NDSolve was exported as an interpolating function, such as:

Export["sol_Tf.wdx",sol]

sol_Tf.wdx is the export file of the NDSolve run.

Also see my previous question.

I then imported it as:

Import["sol_Tf.wdx"]

and used the boundary conditions:

A[x, Ti] == A[x, Tf]
B[x, Ti] == B[x, Tf] 
C[x, Ti] == C[x, Tf]

But the results were different from running it all in one go. What went wrong here? Is there a better way of doing this?

For example, in the heat equation:

sol = NDSolve[{D[u[t, x], t] == D[u[t, x], x, x], u[0, x] == 0,
u[t, 0] == Sin[t], u[t, 5] == 0}, u, {t, 0, 10}, {x, 0, 5}]

I would want to continue it with the end results of the previous run, as the new boundary conditions of the next run; RAM free, and by running it again from t = 10 to t = 20.

Answer 1

You may use the NDSolve Components and Data Structures tutorial to control the memory usage and save down the state of intermediate runs.

Initialise the NDSolve`StateData for the complete range you need to create solutions. Below is done for 0 <= t <= 30 and I will iterate in chunks of 10.

ndsStateData = First@NDSolve`ProcessEquations[
    {
     D[u[t, x], t] == D[u[t, x], x, x],
     u[0, x] == 0,
     u[t, 0] == Sin[t],
     u[t, 5] == 0
     },
    u,
    {t, 0, 30}, {x, 0, 5}
    ];

Next I iterate in 3 chunks of 10. ndsStateData can be saved down after the current chunk's solution is extracted with NDSolve`ProcessSolutions. Below I reassign ndsStateData with the reinitialised NDSolve`StateData instead.

sols = {};
Module[{step = #},
    NDSolve`Iterate[ndsStateData, {(step - 1)*10, step 10}];
    AppendTo[sols, u /. NDSolve`ProcessSolutions[ndsStateData]];
    ndsStateData =
     First@NDSolve`Reinitialize[ndsStateData, {u[step 10, x] == sols[[step]][step 10, x]}];
    ] & /@ Range@3;

sols contains the 3 solutions over each chunk. Notice the difference in the domains.

sols

Plotting the solutions.

Plot3D[sols[[#]][t, x], {t, (# - 1) 10, # 10}, {x, 0, 5},
    PlotRange -> Full,
    PlotStyle -> ColorData[109][#]
    ] & /@ Range@3 // Show[#, PlotRange -> All] &

Hope this helps.

Answer 2

The following throws some inconsistency warnings, but appears to work (i.e., no obvious discontinuities at the seam). Also, I did not change the boundary condition, but I fed the final condition of the first solution to the "initial" condition of the second solution.

uif1 = NDSolveValue[{D[u[t, x], t] == D[u[t, x], x, x], u[0, x] == 0, 
    u[t, 0] == Sin[t], u[t, 5] == 0}, u, {t, 0, 10}, {x, 0, 5}];
uif2 = NDSolveValue[{D[u[t, x], t] == D[u[t, x], x, x], 
    u[10, x] == uif1[10, x], u[t, 0] == Sin[t], u[t, 5] == 0}, 
   u, {t, 10, 20}, {x, 0, 5}];
plt1 = Plot3D[uif1[t, x], {t, 0., 10.}, {x, 0.`, 5.`}, 
   PlotRange -> {{0, 20}, {0, 5}, {-1, 1}}];
plt2 = Plot3D[uif2[t, x], {t, 10.`, 20.`}, {x, 0.`, 5.`}, 
   PlotRange -> {{0, 20}, {0, 5}, {-1, 1}}, 
   PlotStyle -> RGBColor[0.22`, 1.`, 0.44`]];
Show[{plt1, plt2}]