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

Question

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

Answer 1

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 NDSolve`StateData[] 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 = NDSolve`InvokeMethod[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 = NDSolve`SolutionDataComponent[sd, "T"];
   If[IntegerQ[every] && Mod[monstep, every] == 0,    (* save nth step *)
    savevar = NDSolve`ProcessSolutions[state];
    DumpSave[
     file <> If[which === All,
        ToString[monstep],
        ""] <> ".mx" // (ssPrint[#]; #) &,
     savevar],
    If[Round[lasttime/every] < Round[montime/every],  (* save periodically *)
     savevar = NDSolve`ProcessSolutions[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}}  *)