NDSolve`Reinitialize start from previous solution but with different boundary conditions

Question

I am solving a set of coupled partial differential equations (some variation of a 1 dimensional compressible fluid flow). I want to solve the equations (NDSolve, this works fine) but then continue solving with different boundary conditions. I am trying to use NDSolve`Reinitialize for this but it gives me an error.

Current code: I first have someFunction as boundary condition. Then I want to solve the system, and continue with NEWFUNCTION as boundary condition.

(* rhostart, vstart, estart are defined previously *)
boundaryconditions = (
  rho[0, z] == rhostart[z]
 && v[0, z] == vstart[z]
 && e[0, z] == estart[z]
 && ...
 && e[t, 0] == someFunction[t]
 && ...);

statedata = First[NDSolve`ProcessEquations[
    {eq1, eq2, eq3, boundaryconditions},
    {rho, v, e}, {t, 0, tmax}, {z, 0, zmax}
]]

(* I call Iterate in a loop to check if the result is sufficiently stable (I omitted that code here) *)
NDSolve`Iterate[statedata, curTime]

(* When stable, I process the result *)
solution = NDSolve`ProcessSolutions[statedata]
rhostart[z_] := Evaluate[rho[curTime,z] /. solution];
  estart[z_] := Evaluate[e[curTime,z] /. solution];
  vstart[z_] := Evaluate[v[curTime,z] /. solution];

(* I compute a few things using the result, omitted here *)

(* Now I want to change a boundary condition by changing e[t,0] *)
newboundaryconditions = (
  rho[0, z] == rhostart[z]
 && v[0, z] == vstart[z]
 && e[0, z] == estart[z]
 && ...
 && e[t, 0] == NEWFUNCTION[t]
 && ...);
NDSolve`Reinitialize[statedata, {newboundaryconditions}];

This gives me the following error:

NDSolve`Reinitialize::ndsv: Cannot find starting value for the variable v.

More details: One of my boundary conditions involves setting a variable at x=0 to a certain value. I then let the system evolve in time untill it reaches a steady state. Subsequently I want to compute a few values and then change this boundary condition to a different value and let the system evolve again in time starting from where it left off untill it reaches a sufficiently steady state again.

Important: I could change the original boundary condition function someFunction to incorporate all changes at long time intervals. However I do not want to keep the full solution of NDSolve in memory all the time (it already reaches 6 GB quickly) so I think it is neccesary to restart the NDSolve process each time.

If anything is unclear or if more information is needed please let me know.

Answer 1

I can confirm this problem appears when altering the boundary conditions, but not when altering the initial conditions. Here's an example of the former (the latter case is handled successfully in the documentation):

(* setup simple damped wave equation *)
diffeq = Derivative[0, 2][u][x, t] + Derivative[0, 1][u][x, t] - 
   Derivative[2, 0][u][x, t] == 0;
l = 5; tmax = 5; (* set bounds *)
(* set boundary and initial conditions *)
boundary1 = (u[0, t] == 0 && u[l, t] == Sin[t]^2 && u[x, 0] == 0 && 
   Derivative[0, 1][u][x, 0] == 0);
(* create NDSolve data object *)
ndssdata = 
 First@NDSolve`ProcessEquations[{diffeq, boundary1}, 
   u, {x, 0, l}, {t, 0, tmax}];
(* integrate and get solution*)
NDSolve`Iterate[ndssdata, tmax];
soln1 = NDSolve`ProcessSolutions[ndssdata];
(* setup new initial conditions *)
u0 = FunctionInterpolation[u[x, tmax] /. soln1, {x, 0, l, 0.01}][x];
du0 = FunctionInterpolation[
  Derivative[0, 1][u][x, tmax] /. soln1, {x, 0, l, 0.01}][x];
(* set initial conditions and new boundary conditions *)
boundary2 = (u[0, t] == 0 && u[l, t] == Sin[t + tmax]^2 && 
   u[x, 0] == u0 && Derivative[0, 1][u][x, 0] == du0);
(* error on reinitialize *)
ndssdata2 = First@NDSolve`Reinitialize[ndssdata, boundary2]
(* NDSolve`Reinitialize::ndsv: Cannot find starting value for the variable u. *)

As a work-around, you can just re-process the equations. If your individual integrations take a long time (I suspect they do if they generate ~6GB of data), and your equations are not hideously complex, then this should add minimal overhead.

ndssdata2 = 
 First@NDSolve`ProcessEquations[{diffeq, boundary2}, 
   u, {t, 0, tmax}, {x, 0, l}];
NDSolve`Iterate[ndssdata2, tmax];
soln2 = NDSolve`ProcessSolutions[ndssdata2];
(* solutions match! *)
Plot3D[If[t <= tmax, u[x, t] /. soln1, u[x, t - tmax] /. soln2], {x, 
  0, l}, {t, 0, 2 tmax}]

If you're doing it this way, and you don't need the extra control over the solving process, then it's a bit simpler just to use regular NDSolve:

soln1 = First@NDSolve[{diffeq, boundary1}, u, {x, 0, l}, {t, 0, tmax}];
(* set u0, du0, and boundary2 as before *)
soln2 = First@NDSolve[{diffeq, boundary2}, u, {x, 0, l}, {t, 0, tmax}];
solnfull = {u -> 
   FunctionInterpolation[
    If[t <= tmax, u[x, t] /. soln1, u[x, t - tmax] /. soln2], {x, 0, 
     l}, {t, 0, 2 tmax}]};
Plot3D[u[x, t] /. solnfull, {x, 0, l}, {t, 0, 2 tmax}]