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

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} *)
  NDSolve`ProcessEquations @@ (ndscall /.
     {t, a_, b_} :> {t, 
       tinitial = a; b, tfinal = b});
mem1 = MaxMemoryUsed[];
Monitor[
 sol = Table[
    NDSolve`Iterate[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.