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

Question

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?

Answer 1

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.