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/
        kllarge}]\ )\) Cos[\((2\ Pi/
         domainlength)\) RandomInteger[{2, 10}] x] + \((1/
        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"}]]
Table[Flatten[{x, H[x, 30000]} /. solution], {x, 0, domainlength, 


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 /: 
  NDSolve`InitializeMethod[MonitorMethod, stepmode_, sd_, rhs_, 
   state_, OptionsPattern[MonitorMethod]] := Module[{submethod, mf},
   mf = OptionValue["MonitorFunction"];
   submethod = OptionValue[Method];
   If[submethod === Automatic, submethod = "StiffnessSwitching"];
   submethod = 
     NDSolve`InitializeSubmethod[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]];

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


Check the directory:

FileNames["sol*", $TemporaryDirectory]

Retrieve a solution with Get:

<< /private/var/folders/9d/68khy4s15sjf9qfpnhqz9tnc0000gr/T/sol1735.mx

enter image description here

Note that the time domain is around 1735:

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.