# How to improve FDM solver for unsteady viscous flow?

To solve the problem that is discussed in the paper Finite Difference Analysis of Time-Dependent Viscous Nanofluid Flow Between Parallel Plates we developed FDM solver based on the code from the blog Using Mathematica to Simulate and Visualize Fluid Flow in a Box. To check the results obtained, we simultaneously used FEM solver published here. The problem we want to discuss is that FEM solver is 10 times faster than FDM solver. Whereas we expected that there would be an inverse relationship. Since we asked before, why is FEM solver so slow? Let's compare the 2 methods FEM and FDM, and the results. The problem we want to solve is shown in the Figure 1 below

1. FDM solver. So we see that in the equations for flow there is no pressure although the equation of continuity must be satisfied. Therefore we add to the left part of the second and third equation $$\nabla_xP,\nabla_yP$$ respectively with initial condition $$p(0,\xi ,\eta)=0$$ and boundary condition $$p(\tau,l_b,\eta)=0$$, where $$l_b=5$$, and $$0\le \tau \le 1,0 \le \xi \le l_b, 0 \le \eta \le 1$$. Also below we us $$t, x, y$$ instead of $$\tau, \xi, \eta$$. First 3 equations can be solved separately as follows

XYgrid[dom_List, pts_List] :=
N@Range[Sequence @@ #1, Abs[Subtract @@ #1]/#2] &, {dom,
pts - 1}];
BoundaryIndex[xgridlen_, ygridlen_] :=
Module[{tmp, left, right, bot, top},
tmp = Table[(n - 1) ygridlen + Range[1, ygridlen], {n, 1,
xgridlen}]; {left, right} = tmp[[{1, -1}]]; {bot, top} =
Transpose[{First[#], Last[#]} & /@ tmp]; {top, right[[2 ;; -2]],
bot, left[[2 ;; -2]]}];
Attributes[MakeVariables] = {Listable};
MakeVariables[var_, n_] := Table[Unique[var], {n}];
FDMat[deriv_, xygrid_, difforder_] :=
Map[NDSolveFiniteDifferenceDerivative[#, xygrid,
"DifferenceOrder" -> difforder]["DifferentiationMatrix"] &, deriv]
{dt, tmax, domain, pts, difforder, Rey, H1, c1, M1, b1} = {1/20,
20, {{0, 5}, {0, 1}}, {50, 10}, 4, 1., .1, .5, 5, 3};
xygrid = XYgrid[domain, pts]; {nx, ny} =
Map[Length, xygrid]; {top, right, bot, left} =
BoundaryIndex[nx, ny]; {dx, dy, dx2, dy2} =
FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}}, xygrid, difforder]; {dx1,
dy1} = FDMat[{{1, 0}, {0, 1}}, xygrid, 1]; {uvar, vvar, pvar} =
MakeVariables[{u, v, p}, nx ny]; {uvar0, vvar0, pvar0} =
Array[#, nx ny] & /@ {u0, v0, p0}; eqnu = (uvar - uvar0)/dt +
uvar0 (dx . uvar) + vvar0 (dy . uvar) + dx1 . pvar -
1/Rey (dx2 + dy2) .
uvar + (H1 + c1)/(1/M1) uvar; eqnv = (vvar - vvar0)/dt +
uvar0 (dx . vvar) + vvar0 (dy . vvar) + dy1 . pvar -
1/Rey (dx2 + dy2) . vvar + c1/(1/M1) vvar;
eqncont = dx . uvar + dy . vvar; boundaries =
Join[top, right, bot, left]; sgrid =
Flatten[Outer[List, Sequence @@ xygrid], 1];
eqnu[[boundaries]] = uvar[[boundaries]];
eqnu[[bot]] = uvar[[bot]] - M1 sgrid[[bot]][[All, 1]];
eqnu[[right]] = (dx . uvar)[[right]];
eqnv[[boundaries]] = vvar[[boundaries]];
eqnv[[bot]] = vvar[[bot]] - b1; eqnv[[right]] = (dx . vvar)[[right]];
eqncont[[right]] =
pvar[[right]]; {deqns, dvars} = {Join[eqnu, eqnv, eqncont],
Join[uvar, vvar, pvar]};

rule[0] =
Join[Table[uvar0[[i]] -> 0, {i, Length[uvar0]}],
Table[vvar0[[i]] -> 0, {i, Length[vvar0]}]]; dvar0 =
Table[1/10, Length[dvars]];
Do[sol[j] =
Quiet@FindRoot[deqns /. rule[j - 1],
Table[{dvars[[i]], dvar0[[i]]}, {i, Length[dvars]}],
Method -> {"Newton", "StepControl" -> "TrustRegion"}];
dvar0 = dvars /. sol[j];
rule[j] =
Join[Table[uvar0[[i]] -> dvar0[[i]], {i, Length[uvar0]}],
Table[vvar0[[i]] -> dvar0[[i + Length[uvar0]]], {i, 1,
Length[vvar0]}]];, {j, tmax}]; // AbsoluteTiming


On my notebook with 9th gen Intel core 7 it takes 291.367 s. Visualization

{u, v, p} =
Map[Interpolation@Join[sgrid, Transpose@List@#, 2] &,
Partition[dvars /. sol[tmax], Length[sgrid]]];
{ContourPlot[u[x, y], {x, 0, 5}, {y, 0, 1}, ColorFunction -> "Rainbow",
Contours -> 20, PlotLegends -> Automatic, AspectRatio -> Automatic,
MaxRecursion -> 2, PlotPoints -> 50, PlotLabel -> "u",
FrameLabel -> Automatic],

ContourPlot[v[x, y], {x, 0, 5}, {y, 0, 1}, ColorFunction -> "Rainbow",
Contours -> 20, PlotLegends -> Automatic, AspectRatio -> Automatic,
PlotLabel -> "v", FrameLabel -> Automatic]}


2. FEM solver.

UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
P0[0][x_, y_] := 0;
AbsoluteTiming[
Do[{UX[i], VY[i], P0[i]} = NDSolveValue[{{Inactive[Div][{{-\[Mu], 0}, {0, -\[Mu]}} . Inactive[Grad][u[x, y], {x, y}], {x, y}] +
Derivative[1, 0][p][x, y] + ((H1 + c1)/(1/M1))*u[x, y] + (u[x, y] - UX[i - 1][x, y])/dt + UX[i - 1][x, y]*D[u[x, y], x] +
VY[i - 1][x, y]*D[u[x, y], y], Inactive[Div][{{-\[Mu], 0}, {0, -\[Mu]}} . Inactive[Grad][v[x, y], {x, y}], {x, y}] +
Derivative[0, 1][p][x, y] + (c1/(1/M1))*v[x, y] + (v[x, y] - VY[i - 1][x, y])/dt + UX[i - 1][x, y]*D[v[x, y], x] +
VY[i - 1][x, y]*D[v[x, y], y], Derivative[1, 0][u][x, y] + Derivative[0, 1][v][x, y]} == {0, 0, 0} /.
{\[Mu] -> 1, M1 -> 5, H1 -> 0.1, c1 -> 0.5, dt -> 1/20}, {DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 0 && y > 0],
DirichletCondition[{u[x, y] == M1*x, v[x, y] == b1}, y == 0 && x > 0.], DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y == L],
DirichletCondition[{p[x, y] == 0}, x == lb]} /. {L -> 1, lb -> 5, M1 -> 5, b1 -> 3, dt -> 1/20}}, {u, v, p},
Element[{x, y}, Rectangle[{0, 0}, {5, 1}]], Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.0025}}], {i, 1, 20}]; ]


On my notebook it takes 16.4886 s. Visualization

{ContourPlot[UX[20][x, y], {x, 0, 5}, {y, 0, 1},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
Contours -> 20, PlotLabel -> "u", FrameLabel -> Automatic,
AspectRatio -> Automatic, PlotRange -> All, ContourStyle -> Yellow],
ContourPlot[VY[20][x, y], {x, 0, 5}, {y, 0, 1},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
Contours -> 20, PlotLabel -> "v", FrameLabel -> Automatic,
AspectRatio -> Automatic, PlotRange -> All, ContourStyle -> Yellow]}


Therefore, with the same results FEM solver 17-18 times faster than FDM solver. Now we need to make parametric research, for example, compute the Nusselt number $$Nu_x=-x(\frac{\partial \theta}{\partial x}+\frac{\partial \theta}{\partial y})$$ as a function of $$Pr, Nb$$ with a given flow. For this we organize Module as follows

1 FDM solver.

Nu[Pr_, Nb_, x0_, y0_] :=
Module[{Sc1 = 0.5, Nt1 = 0.5}, {dx, dy, dx2, dy2} =
FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}}, xygrid, 2]; {Tvar,
phivar} =
MakeVariables[{\[Theta], \[Phi]}, nx ny]; {Tvar0, phivar0} =
Array[#, nx ny] & /@ {T0, phi0};
eqnT = (Tvar - Tvar0)/dt + uvar0 (dx . Tvar) + vvar0 (dy . Tvar) -
1/Pr (dx2 + dy2) . Tvar -
Nt1 ((dx . Tvar) (dx . Tvar0) + (dy . Tvar) (dy . Tvar0)) -
Nb ((dx . Tvar) (dx . phivar0) + (dy . Tvar) (dy . phivar0));
eqnphi = (phivar - phivar0)/dt + uvar0 (dx . phivar) +
vvar0 (dy . phivar) - 1/Sc1 (dx2 + dy2) . phivar -
Nt1/(Nb*Sc1) (dx2 + dy2) . Tvar;
eqnT[[boundaries]] = Tvar[[boundaries]];
eqnT[[bot]] = Tvar[[bot]] - 1;
eqnT[[right]] = (dx . Tvar)[[right]];
eqnT[[left]] = (dx . Tvar)[[left]];
eqnphi[[boundaries]] = phivar[[boundaries]];
eqnphi[[bot]] = phivar[[bot]] - 1;
eqnphi[[right]] = (dx . phivar)[[right]];
eqnphi[[left]] = (dx . phivar)[[left]]; {deqns1,
dvars1} = {Join[eqnT, eqnphi], Join[Tvar, phivar]};
rule1[0] =
Join[rule[0], Table[Tvar0[[i]] -> 0, {i, Length[Tvar0]}],
Table[phivar0[[i]] -> 0, {i, Length[phivar0]}]];
dvar01 = Table[1/10, Length[dvars1]];
Do[sol1[j] =
Quiet@FindRoot[deqns1 /. rule1[j - 1],
Table[{dvars1[[i]], dvar01[[i]]}, {i, Length[dvars1]}]];
dvar01 = dvars1 /. sol1[j];
rule1[j] =
Join[rule[j],
Table[Tvar0[[i]] -> dvar01[[i]], {i, Length[Tvar0]}],
Table[phivar0[[i]] -> dvar01[[i + Length[Tvar0]]], {i, 1,
Length[phivar0]}]];, {j, tmax}]; {T, phi} =
Map[Interpolation@Join[sgrid, Transpose@List@#, 2] &,
Partition[dvars1 /. sol1[tmax],
Length[sgrid]]]; -x0 (Derivative[1, 0][T][x0, y0] +
Derivative[0, 1][T][x0, y0])];


We test this code in one point

Nu[1, .5, 1, 0] // AbsoluteTiming

Out[]= {117.234, 0.225568}


2 FEM solver.

mesh = UX[20]["ElementMesh"]
Nu[Pr1_, Nb1_, x0_, y0_] :=
Module[{Sc1 = 0.5, Nt1 = 0.5}, T0[0][x_, y_] := 0;
Phi[0][x_, y_] := 0;
Do[
{T0[i], Phi[i]} =
NDSolveValue[{{Inactive[
Div][({{-\[Mu]1, 0}, {0, -\[Mu]1}} .
y}] + (T[x, y] - T0[i - 1][x, y])/dt +
UX[i - 1][x, y]*D[T[x, y], x] +
VY[i - 1][x, y]*
D[T[x, y],
y] - (Nb1*(D[T[x, y], x]*D[Phi[i - 1][x, y], x] +
D[T[x, y], y]*D[Phi[i - 1][x, y], y]) +
Nt1*((D[T0[i - 1][x, y], x])^2 + (D[T0[i - 1][x, y],
y])^2)),
Inactive[
Div][({{-\[Mu]2, 0}, {0, -\[Mu]2}} .
Inactive[Grad][phi[x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-\[Mu]3, 0}, {0, -\[Mu]3}} .
y}] + (phi[x, y] - Phi[i - 1][x, y])/dt +
UX[i - 1][x, y]*D[phi[x, y], x] +
VY[i - 1][x, y]*D[phi[x, y], y]} == {0, 0} /. {\[Mu]1 ->
1/Pr1, \[Mu]2 -> 1/Sc1, \[Mu]3 -> Nt1/(Nb1*Sc1),
dt -> 1/20}, {DirichletCondition[{T[x, y] == 0,
phi[x, y] == 0}, y == L && x > 0],
DirichletCondition[{T[x, y] == 1 , phi[x, y] == 1},
y == 0 && x > 0.]
} /. {L -> 1, lb -> 5}}, {T, phi}, {x, y} \[Element] mesh,
Method -> {"FiniteElement",
"InterpolationOrder" -> {T -> 2, phi -> 2}}];, {i, 1,
20}]; -x0 (Derivative[1, 0][T0[20]][x0, y0] +
Derivative[0, 1][T0[20]][x0, y0])];


Test in one point

Nu[1, .5, 1, 0] // AbsoluteTiming

Out[]= {13.8536, 0.229032}


Therefore, we see that the FEM solver is 8-18 times faster than the FDM solver. And this is very strange, since it seems that FDM should be faster than FEM.

Update 1. We can reduce the integration time by a factor of 6-7 using LeastSquares instead of FinfRoot in FDM code as follows

XYgrid[dom_List, pts_List] :=
N@Range[Sequence @@ #1, Abs[Subtract @@ #1]/#2] &, {dom,
pts - 1}];
BoundaryIndex[xgridlen_, ygridlen_] :=
Module[{tmp, left, right, bot, top},
tmp = Table[(n - 1) ygridlen + Range[1, ygridlen], {n, 1,
xgridlen}]; {left, right} = tmp[[{1, -1}]]; {bot, top} =
Transpose[{First[#], Last[#]} & /@ tmp]; {top, right[[2 ;; -2]],
bot, left[[2 ;; -2]]}];
Attributes[MakeVariables] = {Listable};
MakeVariables[var_, n_] := Table[Unique[var], {n}];
FDMat[deriv_, xygrid_, difforder_] :=
Map[NDSolveFiniteDifferenceDerivative[#, xygrid,
"DifferenceOrder" -> difforder]["DifferentiationMatrix"] &, deriv]
{dt, tmax, domain, pts, difforder, Rey, H1, c1, M1, b1} = {1/20,
20, {{0, 5}, {0, 1}}, {50, 10}, 4, 1., .1, .5, 5, 3};
xygrid = XYgrid[domain, pts]; {nx, ny} =
Map[Length, xygrid]; {top, right, bot, left} =
BoundaryIndex[nx, ny]; {dx, dy, dx2, dy2} =
FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}}, xygrid, difforder]; {dx1,
dy1} = FDMat[{{1, 0}, {0, 1}}, xygrid, 1]; {uvar, vvar, pvar} =
MakeVariables[{u, v, p}, nx ny]; {uvar0, vvar0, pvar0} =
Array[#, nx ny] & /@ {u0, v0, p0}; eqnu = (uvar - uvar0)/dt +
uvar0 (dx . uvar) + vvar0 (dy . uvar) + dx1 . pvar -
1/Rey (dx2 + dy2) .
uvar + (H1 + c1)/(1/M1) uvar; eqnv = (vvar - vvar0)/dt +
uvar0 (dx . vvar) + vvar0 (dy . vvar) + dy1 . pvar -
1/Rey (dx2 + dy2) . vvar + c1/(1/M1) vvar;
eqncont = dx . uvar + dy . vvar; boundaries =
Join[top, right, bot, left]; sgrid =
Flatten[Outer[List, Sequence @@ xygrid], 1];
eqnu[[boundaries]] = uvar[[boundaries]];
eqnu[[bot]] = uvar[[bot]] - M1 sgrid[[bot]][[All, 1]];
eqnu[[right]] = (dx . uvar)[[right]];
eqnv[[boundaries]] = vvar[[boundaries]];
eqnv[[bot]] = vvar[[bot]] - b1; eqnv[[bot[[1]]]] = vvar[[bot[[1]]]];
eqnv[[right]] = (dx . vvar)[[right]];
eqncont[[right]] =
pvar[[right]]; {deqns, dvars} = {Join[eqnu, eqnv, eqncont],
Join[uvar, vvar, pvar]};

rule[0] =
Join[Table[uvar0[[i]] -> 0, {i, Length[uvar0]}],
Table[vvar0[[i]] -> 0, {i, Length[vvar0]}]];
Do[{vec, mat} = CoefficientArrays[deqns /. rule[j - 1], dvars];
sol[j] = LeastSquares[mat, -vec, Method -> "Direct"];
dvar0 = sol[j];
rule[j] =
Join[Table[uvar0[[i]] -> dvar0[[i]], {i, Length[uvar0]}],
Table[vvar0[[i]] -> dvar0[[i + Length[uvar0]]], {i, 1,
Length[vvar0]}]];, {j, tmax}]; // AbsoluteTiming

(*Out[]= {45.9387, Null}*)


Update 2 We can reduce computation time in a case of FEM code using ramp function described in tutorials as follows

rules = {length -> 5, height -> 1, M1 -> 5, b1 -> 3, H1 -> 0.1,
c1 -> .5, \[Mu] -> 1, \[Rho] -> 1};
\[CapitalOmega] = Rectangle[{0, 0}, {length, height}] /. rules;
region = RegionPlot[\[CapitalOmega], AspectRatio -> Automatic]

ClearAll[op, ic, bcs, \[Rho], \[Mu]]
op = {\[Rho]*D[u[t, x, y], t] +
Inactive[
Div][({{-\[Mu], 0}, {0, -\[Mu]}} .
Inactive[Grad][u[t, x, y], {x, y}]), {x,
y}] + \[Rho]*{{u[t, x, y], v[t, x, y]}} .
Inactive[Grad][u[t, x, y], {x, y}] +
D[p[t, x, y], x], \[Rho]*D[v[t, x, y], t] +
Inactive[
Div][({{-\[Mu], 0}, {0, -\[Mu]}} .
Inactive[Grad][v[t, x, y], {x, y}]), {x,
y}] + \[Rho]*{{u[t, x, y], v[t, x, y]}} .
Inactive[Grad][v[t, x, y], {x, y}] + D[p[t, x, y], y],
D[u[t, x, y], x] + D[v[t, x, y], y]}; source = {(c1 + H1)*M1*
u[t, x, y], c1*M1*v[t, x, y], 0};
rampFunction[min_, max_, c_, r_] :=
Function[t, (min*Exp[c*r] + max*Exp[r*t])/(Exp[c*r] + Exp[r*t])]
sf = rampFunction[0, 1, 4, 5]; ic = {u[-1, x, y] == 0,
v[-1, x, y] == 0,
p[-1, x, y] ==
0}; bcs = {DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0},
x == 0],
DirichletCondition[{u[t, x, y] == sf[10 t + 5] M1 x,
v[t, x, y] == sf[10 t + 5] b1}, y == 0 && 0 < x < length],
DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0},
y == height && 0 <= x < length],
DirichletCondition[p[t, x, y] == 0., x == length]} /. rules;
AbsoluteTiming[{xVel, yVel, pressure} =
NDSolveValue[{op + source == {0, 0, 0} /. rules, bcs, ic}, {u, v,
p}, {x, y} \[Element] \[CapitalOmega], {t, -1, 1},
Method -> {"TimeIntegration" -> {"IDA",
"MaxDifferenceOrder" -> 2},
"PDEDiscretization" -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> True,
"SpatialDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}}];]


In this case we have Out[]= {6.50973, Null}. But there is the problem with this approach since computation time for the Nusselt number increases 100 times compare to function Nu[] for FEM introduced above. Really, we have

mesh = xVel["ElementMesh"]

ClearAll[ic1, op1, bc1, T, phi]

ic1 = {T[-1, x, y] == 0,
phi[-1, x, y] ==
0}; op1 = {Inactive[
Div][({{-\[Mu]1, 0}, {0, -\[Mu]1}} .
Inactive[Grad][T[t, x, y], {x, y}]), {x, y}] +
D[T[t, x, y], t] + {{xVel[t, x, y], yVel[t, x, y]}} .
T[t, x, y], {x,
y}] - (Nb1*(Inactive[Grad][T[t, x, y], {x, y}] .
Inactive[Grad][phi[t, x, y], {x, y}]) +
Nt1*(Inactive[Grad][T[t, x, y], {x, y}] .
Inactive[
Div][({{-\[Mu]2, 0}, {0, -\[Mu]2}} .
Inactive[Grad][phi[t, x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-\[Mu]3, 0}, {0, -\[Mu]3}} .
Inactive[Grad][T[t, x, y], {x, y}]), {x, y}] +
D[phi[t, x, y], t] + {{xVel[t, x, y], yVel[t, x, y]}} .
phi[t, x, y], {x, y}]}; bc1 = {DirichletCondition[{T[t, x, y] ==
0, phi[t, x, y] == 0}, y == height && x > 0],
DirichletCondition[{T[t, x, y] == sf[10 t + 5],
phi[t, x, y] == sf[10 t + 5]}, y == 0 && x > 0.]};

Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{T0, Phi} =
NDSolveValue[{Flatten[Activate[op1]] == {0, 0} /. {\[Mu]1 ->
1/Pr1, \[Mu]2 -> 1/Sc1, \[Mu]3 -> Nt1/(Nb1*Sc1)} /.
rules /. {Sc1 -> 0.5, Nt1 -> 0.5, Pr1 -> 1, Nb1 -> .5},
bc1 /. rules, ic1}, {T, phi}, {x, y} \[Element] mesh, {t, -1, 1},
Method -> {"TimeIntegration" -> {"IDA",
"MaxDifferenceOrder" -> 2},
"PDEDiscretization" -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> True,
"SpatialDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {T -> 2, phi -> 2}}}},
EvaluationMonitor :> (currentTime = t;)];]


It takes about 2000s on my laptop, nevertheless it gives the Nusselt number Nu=-(Derivative[0, 1, 0][T0][1, 1, 0] + Derivative[0, 0, 1][T0][1, 1, 0]) of about 0.228415, while with FDM we compute Nu[1, .5, 1, 0]=0.225568, and with linear FEM of about 0.229032.

• @ Alex Trounev: How long does each FindRoot take in your FDM Code? Commented Dec 26, 2023 at 10:09
• @ Alex Trounev: On my machine each iteration takes about 15s to converge. With Method->"ConjugateGradient" this is down to 4s for each iteration. Commented Dec 26, 2023 at 10:47
• @wvt_beginner Thank you very much for your help. With option Method->"ConjugateGradient" we have on every step FindRoot::bdmtd: Value of option Method -> ConjugateGradient is not Automatic, Brent, Secant, AffineCovariantNewton or Newton. Commented Dec 26, 2023 at 12:30
• :I'm on version 13.3 here. What version are you using? Commented Dec 26, 2023 at 12:54
• I also tried your visualisation, but u[x,y] and v[x,y] seem to undefiened in your code. Commented Dec 26, 2023 at 13:08

From the text it's not quite clear if you want to speed up the FEM solver or FDM solver or both. I am going to show you how the FEM solver should be written. The main reason that your FEM solver is slow, is that you manually use an iterative approach for the solver which it can do automatically quite well. Here is how you should write the FEM solver:

Here are the Navier-Stokes equations:

CFDModel = {
y}] + {u[x, y], v[x, y]} . Inactive[Grad][u[x, y], {x, y}] +
Derivative[1, 0][p][x, y],
y}] + {u[x, y], v[x, y]} . Inactive[Grad][v[x, y], {x, y}] +
Derivative[0, 1][p][x, y],
Derivative[0, 1][v][x, y] + Derivative[1, 0][u][x, y]};


source = {+(c1 + H1)*M1*u[x, y], c1*M1*v[x, y], 0};


And now, we solve it, same parameters, same mesh etc:

AbsoluteTiming[
{ufun, vfun, pfun} =
NDSolveValue[{CFDModel + source == {0, 0, 0} /. {\[Mu] -> 1,
M1 -> 5, H1 -> 0.1, c1 -> 0.5,
dt -> 1/20}, {DirichletCondition[{u[x, y] == 0, v[x, y] == 0},
x == 0 && y > 0],
DirichletCondition[{u[x, y] == M1*x, v[x, y] == b1},
y == 0 && x > 0.],
DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y == L],
DirichletCondition[{p[x, y] == 0}, x == lb]} /. {L -> 1,
lb -> 5, M1 -> 5, b1 -> 3, dt -> 1/20}}, {u, v, p},
Element[{x, y}, Rectangle[{0, 0}, {5, 1}]],
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.0025}}];

(* {1.03448, Null} *)


This is more than a factor 10 speedup and this is the documented way for solving the NS-equations.

For completeness:

{ContourPlot[ufun[x, y], {x, 0, 5}, {y, 0, 1},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
Contours -> 20, PlotLabel -> "u", FrameLabel -> Automatic,
AspectRatio -> Automatic, PlotRange -> All, ContourStyle -> Yellow],
ContourPlot[vfun[x, y], {x, 0, 5}, {y, 0, 1},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
Contours -> 20, PlotLabel -> "v", FrameLabel -> Automatic,
AspectRatio -> Automatic, PlotRange -> All, ContourStyle -> Yellow]}


Now, you may argue that in some cases the approach I have shown here does not converge. Sure, that can happen if the nonlinearity is too strong. But even in that case the documented approach is preferable. And the reason is simple: The option "InitialSeeding" can be given to a parameterized version of the Navier-Stokes equations. Since it's not necessary here, I can't quite show it. However, in the upcoming version 14 there is a new tutorial on laminar flow modeling that shows this approach. In the mean time you can have a look at some hyperelastic material models that are very strongly nonlinear and where the same approach is used. Have a look here. You can transfer that approach also to your parametric study.

In summary, simplifying the FEM solver give more that a factor of 10 speedup compared to your FEM solver (and even more compared to your FDM solver).

Update 1:

You have updated your question from a stationary to a transient version. I assume you are satisfied with the performance of the transient Navier-Stokes solver and I thus show the time integration of the second PDE you give:

mesh = xVel["ElementMesh"]

ClearAll[ic1, op1, bc1, T, phi]

ic1 = {T[-1, x, y] == 0, phi[-1, x, y] == 0};
bc1 = {DirichletCondition[{T[t, x, y] == 0, phi[t, x, y] == 0},
y == height && x > 0],
DirichletCondition[{T[t, x, y] == sf[10  t + 5],
phi[t, x, y] == sf[10  t + 5]}, y == 0 && x > 0.]};
op1 = {
Inactive[
Div][({{-\[Mu]1, 0}, {0, -\[Mu]1}} .
Inactive[Grad][T[t, x, y], {x, y}]), {x, y}] +
D[T[t, x, y], t] + {{xVel[t, x, y], yVel[t, x, y]}} .
T[t, x, y], {x,
y}] + ((Nb1*
Grad[T[t, x, y], {x, y}] . Grad[phi[t, x, y], {x, y}]) + (Nt1*
Inactive[
Div][({{-\[Mu]2, 0}, {0, -\[Mu]2}} .
Inactive[Grad][phi[t, x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-\[Mu]3, 0}, {0, -\[Mu]3}} .
Inactive[Grad][T[t, x, y], {x, y}]), {x, y}] +
D[phi[t, x, y], t] + {{xVel[t, x, y], yVel[t, x, y]}} .
Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{T0, Phi} =
NDSolveValue[{op1 == {0, 0} /. {\[Mu]1 -> 1/Pr1, \[Mu]2 ->
1/Sc1, \[Mu]3 -> Nt1/(Nb1*Sc1)} /. rules /. {Sc1 -> 0.5,
Nt1 -> 0.5, Pr1 -> 1, Nb1 -> .5}, bc1 /. rules, ic1}, {T,
phi}, {x, y} \[Element] mesh, {t, -1, 1},
Method -> {"PDEDiscretization" -> {"MethodOfLines"(*,
"DifferentiateBoundaryConditions"->True*)}},
EvaluationMonitor :> (currentTime = t;)];]


For me this time integrates in about 130 seconds. Almost a factor 20 faster than your version. Check if you need the "DifferentiateBoundaryConditions". I have said many times before you never, never, never need to specify "InterpolationOrder" -> {T -> 2, phi -> 2} when the requested interpolation order is the same, this is only useful if at least one of the requested interpolation orders differ. Also, there is no point in Inactivting a PDE and then calling Activate on it.

The one thing that is not ideal here is that I used:

((Nb1*Grad[T[t, x, y], {x, y}] . Grad[phi[t, x, y], {x, y}]) + (Nt1*


in place of the Inactive version. The reason is that the PDE parser stumbles over that; but the solution is not to activate the entire PDE but, rather, just that part. Have a look and see if this solution is what you expect.

• Thank you very much for your approach to speed up FEM solver. But, your solution is valid for steady flow only, while my FEM solver designed for unsteady flow. Also I need to improve FDM solver. Commented Jan 1 at 14:21
• @AlexTrounev same principal for unsteady flow. Commented Jan 1 at 14:23
• Can you show FEM solver for unsteady flow? Commented Jan 1 at 14:38
• @AlexTrounev, you mean like this? Commented Jan 1 at 14:42
• Something like this, but unfortunately this approach is not working. Commented Jan 1 at 15:42

Since my post is too long I would like to answer my question rather then update it. First, I developed new code based on explicit difference scheme with projection step instead of direct pressure computation. Therefore we use on every time step the equations $$\nabla^2P=\nabla.\vec{v}, \vec{v}=\vec{v}-\nabla P$$ instead of continuity equation $$\nabla. \vec{v}=0$$. As result we decrease computation time for FDM code from 291.367 s to 20.6789s that practically is the same as with FEM code (16.4886 s). We have

XYgrid[dom_List, pts_List] :=
N@Range[Sequence @@ #1, Abs[Subtract @@ #1]/#2] &, {dom,
pts - 1}];
BoundaryIndex[xgridlen_, ygridlen_] :=
Module[{tmp, left, right, bot, top, left1, right1},
tmp = Table[(n - 1) ygridlen + Range[1, ygridlen], {n, 1,
xgridlen}]; {left, right} = tmp[[{1, -1}]]; {left1, right1} =
tmp[[{2, -2}]]; {bot, top} =
Transpose[{First[#], Last[#]} & /@ tmp]; {top, right[[2 ;; -2]],
bot, left[[2 ;; -2]], right1[[2 ;; -2]], left1[[2 ;; -2]]}];
Attributes[MakeVariables] = {Listable};
MakeVariables[var_, n_] := Table[Unique[var], {n}];
FDMat[deriv_, xygrid_, difforder_] :=
Map[NDSolveFiniteDifferenceDerivative[#, xygrid,
"DifferenceOrder" -> difforder]["DifferentiationMatrix"] &, deriv]
{dt, tmax, domain, pts, difforder, Rey, H1, c1, M1, b1} = {1/1000,
1000, {{0, 5}, {0, 1}}, {50, 10}, 4, 1., .1, .5, 5, 3};
xygrid = XYgrid[domain, pts]; {nx, ny} =
Map[Length, xygrid]; {top, right, bot, left, right1, left1} =
BoundaryIndex[nx, ny]; {dx, dy, dx2, dy2} =
FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}}, xygrid, difforder]; {dx1,
dy1} = FDMat[{{1, 0}, {0, 1}}, xygrid, 1]; {uvar0, vvar0, pvar0} =
Table[ConstantArray[0, {nx ny, tmax}], {3}]; boundaries =
Join[top, right, bot, left]; sgrid =
Flatten[Outer[List, Sequence @@ xygrid], 1]; {pvar} =
MakeVariables[{pp}, nx ny];

Do[{uvar, vvar} = {uvar0[[All, i]], vvar0[[All, i]]};
eqnu = uvar (dx . uvar) + vvar (dy . uvar) -
1/Rey (dx2 + dy2) . uvar + (H1 + c1)/(1/M1) uvar;
eqnv = uvar (dx . vvar) + vvar (dy . vvar) -
1/Rey (dx2 + dy2) . vvar + c1/(1/M1) vvar;
uvar0[[All, i + 1]] = uvar0[[All, i]] - dt eqnu;
vvar0[[All, i + 1]] = vvar0[[All, i]] - dt eqnv;
uvar0[[top, i + 1]] = uvar0[[top, 1]];
uvar0[[left, i + 1]] = uvar0[[left, 1]];
uvar0[[bot, i + 1]] = M1 sgrid[[bot]][[All, 1]];
vvar0[[bot, i + 1]] = b1; vvar0[[top, i + 1]] = vvar0[[top, 1]];
vvar0[[left, i + 1]] = vvar0[[left, 1]];
eqnp = (dx2 + dy2) . pvar - dx . uvar0[[All, i + 1]] -
dy . vvar0[[All, i + 1]]; eqnp[[right]] = pvar[[right]];
eqnp[[left]] = (dx . pvar)[[left]];
eqnp[[top]] = (dy . pvar)[[top]];
eqnp[[bot]] = (dy . pvar)[[bot]]; {vec, mat} =
CoefficientArrays[eqnp, pvar];
sol = LinearSolve[mat, -vec, Method -> "Krylov"];
uvar0[[All, i + 1]] = uvar0[[All, i + 1]] - dx . sol;
vvar0[[All, i + 1]] = vvar0[[All, i + 1]] - dy . sol;
uvar0[[top, i + 1]] = uvar0[[top, 1]];
uvar0[[left, i + 1]] = uvar0[[left, 1]];
uvar0[[bot, i + 1]] = M1 sgrid[[bot]][[All, 1]];
vvar0[[bot, i + 1]] = b1; vvar0[[top, i + 1]] = vvar0[[top, 1]];
vvar0[[left, i + 1]] = vvar0[[left, 1]], {i, 1,
tmax - 1}]; // AbsoluteTiming


Visualization

{u, v, p} =
Map[Interpolation@Join[sgrid, Transpose@List@#, 2] &,
Partition[Join[uvar0[[All, tmax]], vvar0[[All, tmax]], sol],
Length[sgrid]]];

{ContourPlot[u[x, y], {x, 0, 5}, {y, 0, 1},
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
MaxRecursion -> 2, PlotPoints -> 50, PlotLabel -> "u",
FrameLabel -> Automatic],
ContourPlot[v[x, y], {x, 0, 5}, {y, 0, 1},
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
PlotLabel -> "v"],
ContourPlot[p[x, y], {x, 0, 5}, {y, 0, 1},
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
PlotLabel -> "p"]}


We can compare velocity computed above (dashed lines) with FEM (solid lines) in different sections at t=1

Using velocity we can compute the Nusselt number as follows

Nu[Pr_, Nb_, x0_, y0_] :=
Module[{Sc1 = 0.5, Nt1 = 0.5}, {dx, dy, dx2, dy2} =
FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}}, xygrid, 4]; {Tvar,
phivar} = Table[ConstantArray[0, {nx ny, tmax}], {2}];
Do[eqnT =
uvar0[[All, i]] (dx . Tvar[[All, i]]) +
vvar0[[All, i]] (dy . Tvar[[All, i]]) -
1/Pr (dx2 + dy2) . Tvar[[All, i]] -
Nt1 ((dx . Tvar[[All, i]]) (dx . Tvar[[All, i]]) + (dy .
Tvar[[All, i]]) (dy . Tvar[[All, i]])) -
Nb ((dx . Tvar[[All, i]]) (dx . phivar[[All, i]]) + (dy .
Tvar[[All, i]]) (dy . phivar[[All, i]]));
eqnphi =
uvar0[[All, i]] (dx . phivar[[All, i]]) +
vvar0[[All, i]] (dy . phivar[[All, i]]) -
1/Sc1 (dx2 + dy2) . phivar[[All, i]] -
Nt1/(Nb*Sc1) (dx2 + dy2) . Tvar[[All, i]];

Tvar[[All, i + 1]] = Tvar[[All, i]] - dt eqnT;
phivar[[All, i + 1]] = phivar[[All, i]] - dt eqnphi;
Tvar[[left, i + 1]] = Tvar[[left1, i + 1]];
phivar[[left, i + 1]] = phivar[[left1, i + 1]];
Tvar[[top, i + 1]] = 0; Tvar[[bot, i + 1]] = 1;
phivar[[bot, i + 1]] = 1;
phivar[[top, i + 1]] = 0;, {i, tmax - 1}]; {T, phi} =
Map[Interpolation@Join[sgrid, Transpose@List@#, 2] &,
Partition[Join[Tvar[[All, tmax]], phivar[[All, tmax]]],
Length[sgrid]]]; -x0 (Derivative[1, 0][T][x0, y0] +
Derivative[0, 1][T][x0, y0])];


Test example

Nu[1, .5, 1, 0] // AbsoluteTiming

(*Out[]= {3.17558, 0.22788}*)


Note, that with implicit FDM code we have {117.234, 0.225568}, and with FEM solver we have {13.8536, 0.229032}. Therefore this FDM code for Nu even faster then FEM. Also now we can compile last code for Nu.
Update 1 We can decrease computation time by using mixed explicit-implicit algorithm as follows

XYgrid[dom_List, pts_List] :=
N@Range[Sequence @@ #1, Abs[Subtract @@ #1]/#2] &, {dom,
pts - 1}];
BoundaryIndex[xgridlen_, ygridlen_] :=
Module[{tmp, left, right, bot, top, left1, right1},
tmp = Table[(n - 1) ygridlen + Range[1, ygridlen], {n, 1,
xgridlen}]; {left, right} = tmp[[{1, -1}]]; {left1, right1} =
tmp[[{2, -2}]]; {bot, top} =
Transpose[{First[#], Last[#]} & /@ tmp]; {top, right[[2 ;; -2]],
bot, left[[2 ;; -2]], right1[[2 ;; -2]], left1[[2 ;; -2]]}];
Attributes[MakeVariables] = {Listable};
MakeVariables[var_, n_] := Table[Unique[var], {n}];
FDMat[deriv_, xygrid_, difforder_] :=
Map[NDSolveFiniteDifferenceDerivative[#, xygrid,
"DifferenceOrder" -> difforder]["DifferentiationMatrix"] &, deriv]
{dt, tmax, domain, pts, difforder, Rey, H1, c1, M1, b1} = {1/300,
300, {{0, 5}, {0, 1}}, {50, 10}, 2, 1., .1, .5, 5, 3};
xygrid = XYgrid[domain, pts]; {nx, ny} =
Map[Length, xygrid]; {top, right, bot, left, right1, left1} =
BoundaryIndex[nx, ny]; {dx, dy, dx2, dy2} =
FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}}, xygrid, difforder]; {dx1,
dy1} = FDMat[{{1, 0}, {0, 1}}, xygrid, 1]; {uvar0, vvar0, pvar0} =
Table[ConstantArray[0, {nx ny, tmax}], {3}]; boundaries =
Join[top, right, bot, left]; sgrid =
Flatten[Outer[List, Sequence @@ xygrid], 1]; {pvar} =
MakeVariables[{pp}, nx ny];

matdt1 = IdentityMatrix[nx ny] - dt/Rey (dx2 + dy2); invdt1 =
Inverse[matdt1];

Do[{uvar, vvar} = {uvar0[[All, i]], vvar0[[All, i]]};
eqnu = uvar (dx . uvar) + vvar (dy . uvar) + (H1 + c1)/(1/M1) uvar;
eqnv = uvar (dx . vvar) + vvar (dy . vvar) + (c1)/(1/M1) vvar;
uvar0[[All, i + 1]] = invdt1 . (uvar0[[All, i]] - dt eqnu);
vvar0[[All, i + 1]] = invdt1 . (vvar0[[All, i]] - dt eqnv);
uvar0[[top, i + 1]] = uvar0[[top, 1]];
uvar0[[left, i + 1]] = uvar0[[left, 1]];
uvar0[[bot, i + 1]] = M1 sgrid[[bot]][[All, 1]];
vvar0[[bot, i + 1]] = b1; vvar0[[top, i + 1]] = vvar0[[top, 1]];
vvar0[[left, i + 1]] = vvar0[[left, 1]];
eqnp = (dx2 + dy2) . pvar - dx . uvar0[[All, i + 1]] -
dy . vvar0[[All, i + 1]]; eqnp[[right]] = pvar[[right]];
eqnp[[left]] = (dx . pvar)[[left]];
eqnp[[top]] = (dy . pvar)[[top]];
eqnp[[bot]] = (dy . pvar)[[bot]]; {vec, mat} =
CoefficientArrays[eqnp, pvar];
sol = LinearSolve[mat, -vec, Method -> "Banded"];
uvar0[[All, i + 1]] = uvar0[[All, i + 1]] - dx . sol;
vvar0[[All, i + 1]] = vvar0[[All, i + 1]] - dy . sol;
uvar0[[top, i + 1]] = uvar0[[top, 1]];
uvar0[[left, i + 1]] = uvar0[[left, 1]];
uvar0[[bot, i + 1]] = M1 sgrid[[bot]][[All, 1]];
vvar0[[bot, i + 1]] = b1; vvar0[[top, i + 1]] = vvar0[[top, 1]];
vvar0[[left, i + 1]] = vvar0[[left, 1]], {i, 1,
tmax - 1}]; // AbsoluteTiming


This FDM code takes 3s only, that 100 times less then FDM code with implicit algorithm. This velocity field we compare with FEM solution

{u, v, p} =
Map[Interpolation@Join[sgrid, Transpose@List@#, 2] &,
Partition[Join[uvar0[[All, tmax]], vvar0[[All, tmax]], sol],
Length[sgrid]]];

{ContourPlot[u[x, y], {x, 0, 5}, {y, 0, 1},
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
MaxRecursion -> 2, PlotPoints -> 50, PlotLabel -> "u",
FrameLabel -> Automatic],
ContourPlot[v[x, y], {x, 0, 5}, {y, 0, 1},
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
PlotLabel -> "v"],
ContourPlot[p[x, y], {x, 0, 5}, {y, 0, 1},
ColorFunction -> "Rainbow", Contours -> 20,
PlotLegends -> Automatic, AspectRatio -> Automatic,
PlotLabel -> "p"]}



Also we can compute the Nusselt number as follows

cu0 = Max[Flatten[dx2 + dy2]]/Rey;

Nu[Pr_, Nb_, x0_, y0_] :=
Module[{Sc1 = 0.5, Nt1 = 0.5}, {dx, dy, dx2, dy2} =
FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}}, xygrid, 2];
cu1 = Max[Flatten[dx2 + dy2]];
cu = cu1 Max[{1/Pr, 1/Sc1, Nt1/(Nb*Sc1)}]; dt1 = .7 dt cu0/cu;
tmax1 = Round[tmax dt/dt1]; {Tvar, phivar} =
Table[ConstantArray[0, {nx ny, tmax1}], {2}];
Do[i1 = Round[i dt1/dt];
eqnT = uvar0[[All, i1]] (dx . Tvar[[All, i]]) +
vvar0[[All, i1]] (dy . Tvar[[All, i]]) -
1/Pr (dx2 + dy2) . Tvar[[All, i]] -
Nt1 ((dx . Tvar[[All, i]]) (dx . Tvar[[All, i]]) + (dy .
Tvar[[All, i]]) (dy . Tvar[[All, i]])) -
Nb ((dx . Tvar[[All, i]]) (dx . phivar[[All, i]]) + (dy .
Tvar[[All, i]]) (dy . phivar[[All, i]]));
eqnphi =
uvar0[[All, i1]] (dx . phivar[[All, i]]) +
vvar0[[All, i1]] (dy . phivar[[All, i]]) -
1/Sc1 (dx2 + dy2) . phivar[[All, i]] -
Nt1/(Nb*Sc1) (dx2 + dy2) . Tvar[[All, i]];
Tvar[[All, i + 1]] = Tvar[[All, i]] - dt1 eqnT;
phivar[[All, i + 1]] = phivar[[All, i]] - dt1 eqnphi;
Tvar[[left, i + 1]] = Tvar[[left1, i + 1]];
phivar[[left, i + 1]] = phivar[[left1, i + 1]];
Tvar[[top, i + 1]] = 0; Tvar[[bot, i + 1]] = 1;
phivar[[bot, i + 1]] = 1;
phivar[[top, i + 1]] = 0;, {i, tmax1 - 1}]; {T, phi} =
Map[Interpolation@Join[sgrid, Transpose@List@#, 2] &,
Partition[Join[Tvar[[All, tmax1]], phivar[[All, tmax1]]],
Length[sgrid]]]; -x0 (Derivative[1, 0][T][x0, y0] +
Derivative[0, 1][T][x0, y0])];


Test example

Nu[1, .5, 1, 0] // AbsoluteTiming

Out[]= {2.32057, 0.221814}

• @ Alex Trounev: Nice to see that you figured it out. The pressure field also looks a lot more reasonable. On my machine I get an additional time advantage of about 5.5swith the Method->"Banded" in LinearSolve. Can you confirm this? Commented Jan 6 at 11:15
• @wvt_beginner Yes, you are right. Time decreases from 20.6789s to 18.6551s with Method->"Banded". Maybe we also should add Compile, since LinearSolve is compliable function. Commented Jan 6 at 11:34
• Well, to my knowledge LinearSolve is already highly optimized. By giving it a try with Compile, i get a slowdown of approx. factor 5, with CompilationTarget->"C". Maybe the biggest part of the slowdown is the fact, that - to my knowledge - one cannot use SparseArray objects inside Compile. When using LinearSolve[Normal@mat,(-1)*Normal@vec] I get comparable timings here. Commented Jan 6 at 14:46
• If compilation does not help, then the obvious way is to increase the integration time step. In the implicit scheme we used dt=1/20, and in the code above dt=1/1000. This explicit scheme stable up to dt= 1/500. But code for Nu is unstably in this case. Commented Jan 7 at 1:21