# NDSolve for Fluid Flow - Monitor Residuals

I'm using NDSolve for fluid flow and would like to monitor the convergence of the solution more closely. I have a simple stationary case set up. NDSolve should return the solutions for u,v and p. During the calculation, I want NDSolve to print the current residuals for u, v and p so I can see if the solution actually converges. My first idea was to use StepMonitor, but apparently Mathematica treats non-transient computations as a single step.

Any help on this is greatly appreciated.

(*Fluid*)
\[Rho] = 1;
\[Mu] = 10^-3;

(*Region*)
meshReg =
DiscretizeRegion[
RegionDifference[
Rectangle[{0,0},{2.2,0.41}],
Disk[{1/5,1/5}, 1/20]
],
MaxCellMeasure -> 0.0001
];

(*Operator*)
op = {
D[u[x,y],x]+D[v[x,y],y]
};

(*BoundaryConditions*)
bcs = {
DirichletCondition[{u[x,y] == 0.15, v[x,y] == 0.}, x == 0],
DirichletCondition[{u[x,y] == 0., v[x,y] == 0.}, 0 < x < 2.2],
DirichletCondition[p[x,y] == 0, x == 2.2]
};

(*Solution*)
{time,{xVel, yVel, pressure}} =
AbsoluteTiming@NDSolveValue[
{op == {0,0,0}, bcs}, {u,v,p}, {x,y} \[Element] meshReg,
Method -> {"FiniteElement","InterpolationOrder"->{u->2,v->2,p->1}}
]


I assume you have seen the documentation on the NDSolve options for FEM? There is a section FindRootOptions on how to set up the StepMonitor, EvaluationMonitor and the EvaluationMonitor for the Jacobian evaluation.

You can also see the actual printout of the convergence of the algorithm; this was used while implementing the algorithm and comparing to the reference implementation. For the exact details of what these numbers mean I refer you to references mentioned in the documentation of the Affine corvariant Newton method or in this question.

One last thing before we go: If NDSolve does not complain you can assume that the algorithm has converged to the prescribed tolerance. That is quite different from some other (FEM) tools where it is the customers responsibility to check the convergence. No need to check this here. You will get a message if it does not converge. Don't waste your time on checking that.

NDSolveFEMNonlinearSolversDumpDebugPrintQ[level_] :=
If[TrueQ[level < 2], True, False]

(*Solution*)
{time, {xVel, yVel, pressure}} =
AbsoluteTiming@
NDSolveValue[{op == {0, 0, 0}, bcs}, {u, v, p}, {x, y} \[Element]
meshReg,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]

k=0    ||DX||=0.110323  ||FX||=0.0000397922 \[Lambda]=0.01

k=0 *  ||DX||=0.10922   ||FX||=0.0000393943 \[Lambda]=0.01

k=0 *  ||DX||=0.0233193 ||FX||=2.86501*10^-6    \[Lambda]=1.

k=1    ||DX||=0.0232118 ||FX||=2.86501*10^-6    \[Lambda]=Broyden

k=2    ||DX||=0.015448  ||FX||=2.86501*10^-6    \[Lambda]=1.

k=2 *  ||DX||=0.00198797    ||FX||=3.16805*10^-7    \[Lambda]=1.

k=3    ||DX||=0.00205934    ||FX||=3.16805*10^-7    \[Lambda]=Broyden

k=4    ||DX||=0.000515059   ||FX||=7.05851*10^-8    \[Lambda]=Broyden

k=5    ||DX||=0.000125937   ||FX||=1.62026*10^-8    \[Lambda]=Broyden

k=6    ||DX||=0.0000444832  ||FX||=5.39815*10^-9    \[Lambda]=Broyden

k=7    ||DX||=7.11939*10^-6 ||FX||=1.55357*10^-9    \[Lambda]=Broyden

k=8    ||DX||=2.40416*10^-6 ||FX||=4.43182*10^-10   \[Lambda]=Broyden

k=9    ||DX||=4.11268*10^-7 ||FX||=7.92958*10^-11   \[Lambda]=Broyden

k=10    ||DX||=1.06446*10^-7    ||FX||=2.46871*10^-11   \[Lambda]=Broyden

k=11    ||DX||=2.63202*10^-8    ||FX||=7.16103*10^-12   \[Lambda]=Broyden

k=12    ||DX||=7.71644*10^-9    ||FX||=1.90203*10^-12   \[Lambda]=Broyden
`