How to improve FDM solver for unsteady viscous flow?

Question

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] := 
  MapThread[
   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[NDSolve`FiniteDifferenceDerivative[#, 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}} . 
               Inactive[Grad][T[x, y], {x, y}]), {x, 
              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}} . 
               Inactive[Grad][T[x, y], {x, y}]), {x, 
              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] := 
  MapThread[
   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[NDSolve`FiniteDifferenceDerivative[#, 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], c1M1v[t, x, y], 0};
rampFunction[min_, max_, c_, r_] := 
 Function[t, (minExp[cr] + maxExp[rt])/(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]}} . 
    Inactive[Grad][
     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[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]}} . 
    Inactive[Grad][
     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.

Answer 1

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 = {
   Inactive[Div][-(\[Mu]*Inactive[Grad][u[x, y], {x, y}]), {x, 
      y}] + {u[x, y], v[x, y]} . Inactive[Grad][u[x, y], {x, y}] + 
    Derivative[1, 0][p][x, y], 
   Inactive[Div][-(\[Mu]*Inactive[Grad][v[x, y], {x, y}]), {x, 
      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]};

Here is your source term:

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]}} . 
     Inactive[Grad][
      T[t, x, y], {x, 
       y}] + ((Nb1*
        Grad[T[t, x, y], {x, y}] . Grad[phi[t, x, y], {x, y}]) + (Nt1*
        Grad[T[t, x, y], {x, y}] . 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]}} . 
     Inactive[Grad][phi[t, x, y], {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/(Nb1Sc1)} /. 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*
    Grad[T[t, x, y], {x, y}] . Grad[T[t, x, y], {x, y}]))

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.

Answer 2

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] := 
  MapThread[
   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[NDSolve`FiniteDifferenceDerivative[#, 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] := 
  MapThread[
   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[NDSolve`FiniteDifferenceDerivative[#, 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}