# How to prevent NDSolve to store result for each time step?

It's very painrful especially using FEM that Mathematica gives a solution for transient Navier-Stokes equation at each time step. Is it possicle to control output? This is my code:

Ω = RegionDifference[Rectangle[{-1, -1}, {1, 1}],
RegionUnion[
Flatten[Table[
Disk[{n, m}, 1/16], {n, -1/2, 1/2, 1/2}, {m, -1/2, 1/2, 1/2}]]]];
op = {Derivative[1, 0, 0][u][t, x, y] +
10^-3 Inactive[Div][-Inactive[Grad][u[t, x, y], {x, y}], {x,
y}] + {u[t, x, y], v[t, x, y]}.Inactive[Grad][
u[t, x, y], {x, y}] + Derivative[0, 1, 0][p][t, x, y] - y Tanh[t],
Derivative[1, 0, 0][v][t, x, y] +
10^-3 Inactive[Div][-Inactive[Grad][v[t, x, y], {x, y}], {x,
y}] + {u[t, x, y], v[t, x, y]}.Inactive[Grad][
v[t, x, y], {x, y}] + Derivative[0, 0, 1][p][t, x, y] + x Tanh[t],
Derivative[0, 1, 0][u][t, x, y] + Derivative[0, 0, 1][v][t, x, y]};
bcs =
{DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, True],
DirichletCondition[p[t, x, y] == 0, x == 1 && y == 1]};
ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};
Monitor[
AbsoluteTiming[
{xVel, yVel, pressure} =
NDSolveValue[
{op == {0, 0, 0}, bcs, ic}, {u, v, p}, {x, y} ∈ Ω, {t, 0, 5},
Method ->
{"PDEDiscretization" ->
{"MethodOfLines",
"SpatialDiscretization" ->
{"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.002, "MaxBoundaryCellMeasure" -> 0.02},
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}},
EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];],
currentTime]

If we look at data grid of the solution

Dimensions[xVel[[4]]]
{494, 2, 10492}

So it returns the whole data (values at all mesh points for each time step). It takes a lot of memory if simulation is long. I don't want that Mathematica accomelates results at each time step. I need snapshosts of the flow, say at 10 ( instead of 494) time points. What should I do for that?

• Can you show an example with code? Mar 22, 2021 at 15:51
• Try NDSolve[{x''[t] == x[t]^2, x[0] == 1, x'[0] == 0}, x, {t, 2, 2}]. (Note the time interval.) It's in the docs somewhere, but I can't find it right now. You can also accumulate solution data for a subinterval {t, 1, 2}. Mar 22, 2021 at 15:51
• Apparently I couldn't find it a couple of months ago either: mathematica.stackexchange.com/a/238226/4999 Mar 22, 2021 at 15:57
• NDSolve[{x''[t] == x[t]^2, x[0] == 1, x'[0] == 0}, x[2], {t, 0,2}] returns only the last timestep! Mar 22, 2021 at 17:22
• @user21 My interpretation led me to do the following: Use NDSolveIterate to advance the time integration to a sequence of specified time stops; at each stop, harvest the "SolutionData". This can be used to construct a sequence of interpolations over the spatial discretization, one for each time stop. I wasn't going to worry about interpolation over time unless specifically requested. The sequence would be useful for illustrating the evolution of the system while saving memory. Mar 23, 2021 at 14:26

Sorry just to drop code, but I'm bit too busy at the moment to write up a presentation. Hopefully the code with still be helpful. Ask questions and I'll try to address them later.

The testing setup (testing memory usage). Execute this each time after quitting kernel, before running memory usage tests below:

Ω =
RegionDifference[Rectangle[{-1, -1}, {1, 1}],
RegionUnion[
Flatten[Table[
Disk[{n, m}, 1/16], {n, -1/2, 1/2, 1/2}, {m, -1/2, 1/2, 1/2}]]]];

op = {Derivative[1, 0, 0][u][t, x, y] +
10^-3 Inactive[Div][-Inactive[Grad][u[t, x, y], {x, y}], {x,
y}] + {u[t, x, y], v[t, x, y]} .
Inactive[Grad][u[t, x, y], {x, y}] +
Derivative[0, 1, 0][p][t, x, y] - y Tanh[t],
Derivative[1, 0, 0][v][t, x, y] +
10^-3 Inactive[Div][-Inactive[Grad][v[t, x, y], {x, y}], {x,
y}] + {u[t, x, y], v[t, x, y]} .
Inactive[Grad][v[t, x, y], {x, y}] +
Derivative[0, 0, 1][p][t, x, y] + x Tanh[t],
Derivative[0, 1, 0][u][t, x, y] +
Derivative[0, 0, 1][v][t, x, y]};
bcs = {DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, True],
DirichletCondition[p[t, x, y] == 0, x == 1 && y == 1]};
ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};

(* Hold the NDSolveValue call; release in the two test phases *)
ndscall = Hold[
NDSolveValue
][{op == {0, 0, 0}, bcs, ic}, {u, v,
p}, {x, y} ∈(*emesh*)Ω, {t, 0, 5},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.002,
"MaxBoundaryCellMeasure" -> 0.02},
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}},
EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];

Memory usage of accumulating the whole solution:

Quit[]

(* initialize before executing this cell *)
mem0 = MaxMemoryUsed[];
Monitor[AbsoluteTiming[{xVel, yVel, pressure} =
ReleaseHold[ndscall];], currentTime] (* takes ~100 sec *)
MaxMemoryUsed[] - mem0

(*  560095104  > 183063416 below *)

Memory usage of saving the solution data at t = 0, 1, 2,...5:

Quit[]

(* initialize before executing this cell *)
mem0 = MaxMemoryUsed[];
{state} = (* change time interval from {t, a, b} to {t, b, b} *)
NDSolveProcessEquations @@ (ndscall /.
{t, a_, b_} :> {t,
tinitial = a; b, tfinal = b});
mem1 = MaxMemoryUsed[];
Monitor[
sol = Table[
NDSolveIterate[state, t0];
memused[t0] = MaxMemoryUsed[] - mem1;
state@"SolutionData"["Forward"],
{t0, Subdivide[tinitial, tfinal, 5]}(* time stops *)
];,
currentTime]
MaxMemoryUsed[] - {mem0, mem1}

(*  {183063416, 82331096}  < 560095104 above *)

AssociationMap[memused, Subdivide[tinitial, tfinal, 5]]

(* memory usage once NDSolveIterate loop starts:
Once Iterate starts, each step adds only a few hundred KB
<|0 -> 0,        1 -> 80771560, 2 -> 80771560,
3 -> 81510592, 4 -> 81717312, 5 -> 82331096|>
*)

Processing the solution data. I'm not sure yet that I've figured out the best workflow here. This code processes the last time stop (sol[[-1]]) only, to compare with the final solution

Needs@"NDSolveFEM";
meshes = Through[
NDSolveSolutionDataComponent[state@"Variables", "X"][
"ElementMesh"] /. NDSolveProcessSolutions[state]];
pts = Length /@ Through[meshes@"Coordinates"];
puvVals = TakeList[
NDSolveSolutionDataComponent[sol[[-1]], "X"],
pts];
puvSol = MapThread[#1 -> ElementMeshInterpolation[#2, #3] &,
{NDSolveSolutionDataComponent[state@"Variables", "X"],
meshes,
puvVals}
]

Check last solution data (puvSol) with solution (from ProcessEquations[]):

soltf = NDSolveProcessSolutions[state];
Differences@{v["ValuesOnGrid"] /. soltf // First,
v["ValuesOnGrid"] /. puvSol} // MinMax

(*  {0., 0.}  *)

Summary. Set up the computation with NDSolveProcessEquations and a time integration interval of the form {t, tf, tf}, where tf is the final time. Advance the integration with NDSolveIterate to desired time stops at which the solution is desired. Harvest that data from the state object returned by ProcessEquations.

{state} = NDSolveProcessEquations[{op == {0, 0, 0}, bcs, ic},
{u, v, p}, {x, y} ∈ Ω,
{t, 5, 5},    (***  N.B. No intermediate data saved  ***)
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.002,
"MaxBoundaryCellMeasure" -> 0.02},
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}},
EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])]

sol = Table[
NDSolve`Iterate[state, t0];       (* advance the time integration *)
state@"SolutionData"["Forward"],  (* harvest solution data *)
{t0, {t1, t2, ...}}               (* time stops *)
]

The solution data contains the function values in a flat list, which needs to be subdivided and mapped onto the ElementMesh for the spatial domain.

• This is working sheme. Thanks a lot. Mar 24, 2021 at 21:57