2
$\begingroup$

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?

$\endgroup$
18
  • 1
    $\begingroup$ Can you show an example with code? $\endgroup$
    – MarcoB
    Mar 22, 2021 at 15:51
  • 3
    $\begingroup$ 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}. $\endgroup$
    – Michael E2
    Mar 22, 2021 at 15:51
  • $\begingroup$ Apparently I couldn't find it a couple of months ago either: mathematica.stackexchange.com/a/238226/4999 $\endgroup$
    – Michael E2
    Mar 22, 2021 at 15:57
  • 1
    $\begingroup$ NDSolve[{x''[t] == x[t]^2, x[0] == 1, x'[0] == 0}, x[2], {t, 0,2}] returns only the last timestep! $\endgroup$ Mar 22, 2021 at 17:22
  • 1
    $\begingroup$ @user21 My interpretation led me to do the following: Use NDSolve`Iterate 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. $\endgroup$
    – Michael E2
    Mar 23, 2021 at 14:26

1 Answer 1

4
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ This is working sheme. Thanks a lot. $\endgroup$ Mar 24, 2021 at 21:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.