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 NDSolve`Iterate 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@"NDSolve`FEM`";
meshes = Through[
NDSolve`SolutionDataComponent[state@"Variables", "X"][
"ElementMesh"] /. NDSolve`ProcessSolutions[state]];
pts = Length /@ Through[meshes@"Coordinates"];
puvVals = TakeList[
NDSolve`SolutionDataComponent[sol[[-1]], "X"],
pts];
puvSol = MapThread[#1 -> ElementMeshInterpolation[#2, #3] &,
{NDSolve`SolutionDataComponent[state@"Variables", "X"],
meshes,
puvVals}
]
Check last solution data (puvSol) with solution (from ProcessEquations[]):
soltf = NDSolve`ProcessSolutions[state];
Differences@{v["ValuesOnGrid"] /. soltf // First,
v["ValuesOnGrid"] /. puvSol} // MinMax
(* {0., 0.} *)
Summary. Set up the computation with NDSolve`ProcessEquations and a time integration interval of the form {t, tf, tf}, where tf is the final time. Advance the integration with NDSolve`Iterate to desired time stops at which the solution is desired. Harvest that data from the state object returned by ProcessEquations.
{state} = NDSolve`ProcessEquations[{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.
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}. – Michael E2 Mar 22 '21 at 15:51NDSolveValue[{x''[t] == x[t]^2, x[0] == 1, x'[0] == 0}, {x[1], x[2]}, {t, 0, 2}]gives the solution at two times, accuracy and resolution isn't affected! – Ulrich Neumann Mar 22 '21 at 19:15NDSolve\ProcessEquations(your comment&link) but I don't know how to split the coodinatest` and the spatial coordinates... – Ulrich Neumann Mar 22 '21 at 21:17{t, 0, 100}, there are only 32 steps. The start-up overhead tends to dominate the marginal increase in memory per step. With so few steps, it's hard to decrease the number of steps saved in order to show any significant memory savings. – Michael E2 Mar 23 '21 at 04:06NDSolve`Iterateto 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. – Michael E2 Mar 23 '21 at 14:26ElementMeshInterpolationto get a single interpolating function. But I am not sure about the overall accuracy of this approach, because the time stepper does not store unnecessary data, I would think. But I am interested to see what comes out of this. – user21 Mar 23 '21 at 15:00