# Output of intermediate results from very long NDSolve computation into sequenced CSV files

There is a problem that I was not able to solve for months, if not years, and I can't postpone it any longer. It is this: I know how to output the result of a computation from NDSolve into a file after NDSolve has finished completely, that is, the integration is completed from tinitial to tfinal. But I need to output at some control points of my choosing between tinitial and tfinal, because computation may run for hours and if smth. happens I will loose everything. If I saved a few intermediate results though, then I can re-start from the last one I saved. Here is the simplified code, which I trimmed to run for < 20 sec. I want to insert the output statements similar to the last Export statement in this code into WhenEvent[..] statement that I already have inside NDSolve (presently, it just outputs the t value reached by the integrator) , to dump the results at t=tinitial+deltat, t=tinitial+2deltat, etc. into the files with names like test1.csv, test2.csv, etc. In this example there would be 2 such files, corresponding to t=10000 and 20000. Please help. I am convinced that the solution will help many people. There is nothing remotely close in the posted questions and answers.

RHSH = D[H[x, t] D[H[x, t], x, x], x, x] +
1/2 D[D[H[x, t], x]^2, x, x] +
Subscript[Z, 4] D[H[x, t]^3 D[H[x, t], x, x, x, x], x, x] +
Subscript[Z, 13]
D[H[x, t]^2 D[H[x, t], x] D[H[x, t], x, x, x], x, x];
Nkmax = 0.0947; N\[Lambda]max = 2 Pi/Nkmax; domainlength =
20 N\[Lambda]max; NH0 = 11.47;
klsmall = 10; kllarge = 100; kf = 12; u = 5/2; \[Phi] = 1/2;
IniShapeH = NH0 + \[Phi] Cos[30 (2 Pi/domainlength) x] + u*\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$k = 1$$, $$kf/ 2$$]$$(\((1/ RandomInteger[{klsmall, kllarge}]\ )$$ Cos[$$(2\ Pi/ domainlength)$$ RandomInteger[{2, 10}] x] + $$(1/ RandomInteger[{klsmall, kllarge}]\ )$$ Sin[$$(2\ Pi/ domainlength)$$ RandomInteger[{2, 10}] x])\)\);
\[CapitalDelta]t = 10000; endtime = 30000;
Timing[solution =
NDSolve[{D[H[x, t], t] == RHSH /. {Subscript[Z, 4] -> 0.6322,
Subscript[Z, 13] -> 3.7932}, H[0, t] == H[domainlength, t],
H[x, 0] == IniShapeH,
WhenEvent[Mod[t, \[CapitalDelta]t] == 0, Print[t]]},
H, {x, 0, domainlength}, {t, 0, endtime},
AccuracyGoal -> MachinePrecision/2,
PrecisionGoal -> MachinePrecision/2, InterpolationOrder -> All,
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> 4, MinPoints -> 2500, MaxPoints -> 2500,
AccuracyGoal -> MachinePrecision/2,
PrecisionGoal -> MachinePrecision/2}, Method -> "BDF"}]]
Export["c:\\Test\\test.csv",
Table[Flatten[{x, H[x, 30000]} /. solution], {x, 0, domainlength,
domainlength/1000}]]


## Introduction

Here is a way to save a solution to disk at some point in the middle of an NDSolve integration. It is based on the MonitorMethod in the documentation (see below). While I called the monitor function saveSteps, it writes the current solution (interpolating functions). It uses DumpSave, so you have to give saveSteps a symbol to use. One could use LocalObject instead, if preferred.

## How to save intermediate steps of a time integration (ODE/PDE)

Here is the MonitorMethod from this tutorial, tweaked to pass the NDSolveStateData[] object to the monitor function:

Options[MonitorMethod] = {Method -> Automatic,
"MonitorFunction" -> Function[{h, sd, state, mord},
Print[{"H" -> h, "SD" -> sd, "DifferenceOrder" -> mord}]]};

MonitorMethod /:
NDSolveInitializeMethod[MonitorMethod, stepmode_, sd_, rhs_,
state_, OptionsPattern[MonitorMethod]] := Module[{submethod, mf},
mf = OptionValue["MonitorFunction"];
submethod = OptionValue[Method];
If[submethod === Automatic, submethod = "StiffnessSwitching"];
submethod =
NDSolveInitializeSubmethod[MonitorMethod, submethod, stepmode,
sd, rhs, state];
MonitorMethod[submethod, mf]];

MonitorMethod[submethod_, mf_]["Step"[f_, h_, sd_, state_]] :=
Module[{res = NDSolveInvokeMethod[submethod, f, h, sd, state]},
mf[h, sd, state, submethod["DifferenceOrder"]];
If[SameQ[res[[-1]], submethod], res[[-1]] = Null,
res[[-1]] = MonitorMethod[res[[-1]], mf]];
res]

MonitorMethod[___]["StepInput"] = {"Function"[All], "H",
"SolutionData", "StateData"};
MonitorMethod[___]["StepOutput"] = {"H", "SD", "MethodObject"};

MonitorMethod[submethod_, ___][prop_] := submethod[prop];


Here is a monitor function for saving the current solution at intermediate steps. There are several alternatives. If every is an integer $$n$$, the every $$n$$-th step will be saved. If every is a real number $$\Delta t$$, then the solution will be saved when the integration time step crosses a Mod[t, Δt] == 0 threshold. If which === All, then the filenames will be appended with the step number or a rounded time stamp (with the decimal place omitted in cases where $$\Delta t < 1$$). Caveat: There are some global variables being used. These should be localized in a package, but I felt this unnecessary in a proof-of-concept demonstration.

(*
saveSteps[file, n_Integer, which] save every nth step
saveSteps[file, Δt_Real, which] save step every Δt
saveSteps[file, every, All] save every step in unique file, otherwise last step
*)
ClearAll[saveSteps, iSaveSteps, ssPrint];
ssPrint = Null &; (* ssPrint = Print for verbose output *)
saveSteps[file_String, savevar_Symbol, every_: 1, which_: Last] /;
MatchQ[N@every, _Real] := (
monstep = 0;
montime = 0;
iSaveSteps[StringTrim[file, ".mx"], savevar,
If[IntegerQ[every], every, N@every], which]);
iSaveSteps[file_, savevar_, every_, which_][h_, sd_, state_, dord_] :=
Block[{savevar, lasttime = montime},   (* FAILS if savevar is not a symbol *)
++monstep;
montime = NDSolveSolutionDataComponent[sd, "T"];
If[IntegerQ[every] && Mod[monstep, every] == 0,    (* save nth step *)
savevar = NDSolveProcessSolutions[state];
DumpSave[
file <> If[which === All,
ToString[monstep],
""] <> ".mx" // (ssPrint[#]; #) &,
savevar],
If[Round[lasttime/every] < Round[montime/every],  (* save periodically *)
savevar = NDSolveProcessSolutions[state];
DumpSave[
file <> If[which === All,
If[every < 1,
ToString[Round[montime/10^Floor[Log10[every]]]],
ToString[Round@montime]],
""] <> ".mx" // (ssPrint[#]; #) &, savevar]
]]
];


## Applied to OP's example

Here is the method applied to a simplified version of the OP's problem. I reduced the precision and accuracy goals and the size of the grid, so that it can be tested in less than a second. One significant issue here is that NDSolve tells me that the "BDF" method cannot be used as a submethod.

PrintTemporary@Dynamic@{Clock[Infinity], tmp};  (* progress monitor *)
Block[{ssPrint = Print},  (* verbose output prints filenames as they are created *)
Clear[tempsol];          (* use tempsol to store solutions: must be cleared *)
NDSolve[{
D[H[x, t], t] == RHSH /. {Subscript[Z, 4] -> 0.6322,
Subscript[Z, 13] -> 3.7932}, H[0, t] == H[domainlength, t],
H[x, 0] == IniShapeH},
H, {x, 0, domainlength}, {t, 0, endtime/10},
StepMonitor :> (tmp = t),  (* for progress monitor *)
AccuracyGoal -> MachinePrecision/2,
PrecisionGoal -> MachinePrecision/2, InterpolationOrder -> All,
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> 4, "MinPoints" -> 250, "MaxPoints" -> 250,
AccuracyGoal -> MachinePrecision/4,
PrecisionGoal -> MachinePrecision/4},
Method -> {MonitorMethod,
"MonitorFunction" ->
saveSteps[FileNameJoin[{$TemporaryDirectory, "sol"}], tempsol, 10000./10, All]}}]; ]  /private/var/folders/9d/68khy4s15sjf9qfpnhqz9tnc0000gr/T/sol503.mx /private/var/folders/9d/68khy4s15sjf9qfpnhqz9tnc0000gr/T/sol1735.mx /private/var/folders/9d/68khy4s15sjf9qfpnhqz9tnc0000gr/T/sol2635.mx Check the directory: FileNames["sol*",$TemporaryDirectory]
(*
{"/private/var/folders/9d/68khy4s15sjf9qfpnhqz9tnc0000gr/T/sol1735.mx",
"/private/var/folders/9d/68khy4s15sjf9qfpnhqz9tnc0000gr/T/sol2635.mx",
"/private/var/folders/9d/68khy4s15sjf9qfpnhqz9tnc0000gr/T/sol503.mx"}
*)


Retrieve a solution with Get:

<< /private/var/folders/9d/68khy4s15sjf9qfpnhqz9tnc0000gr/T/sol1735.mx
tempsol Note that the time domain is around 1735:

H["Domain"] /. tempsol
(*  {{0., 1326.97}, {0., 1734.91}}  *)
`