Speeding up Boussinesq equations solving?

Question

I am working on Boussinesq equation. The notebook can run perfectly for only 0.6 steps and then the calculation starts running slowly after 0.7. All boundary conditions seemed fine. I am unsure if I used the correct method to solve this problem.

   tmax = 10;
Monitor[
 AbsoluteTiming[
  {xVel, yVel, pressure, temperature} = NDSolveValue[
     Flatten@{op == {0, 0, 0, 
         Subscript[\[CapitalGamma], convectiveLeft] + 
          Subscript[\[CapitalGamma], convectiveRight] + 
          Subscript[\[CapitalGamma], convectiveTop]}, bcsflow, 
       bcstemperatures, ic}, {u, v, p, T}, {x, y} \[Element] mesh, {t,
       0, tmax},
     Method -> {"TimeIntegration" -> {"IDA", 
         "MaxDifferenceOrder" -> 2, 
         "ImplicitSolver" -> {"Newton", 
           "LinearSolveMethod" -> "Pardiso"}},
       "PDEDiscretization" -> {
         "MethodOfLines", "SpatialDiscretization" -> {
           "FiniteElement",
           "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}}
         }
       },
     EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])
     ];
  ],
 currentTime
 ]

These are my error msg.

NDSolveValue::fembcib: The finite element mesh has internal boundaries. Boundary conditions with 'True' as a predicate will set boundary conditions at these internal boundaries. It is recommended to use a less general predicate than 'True' for the boundary conditions {DirichletCondition[p==0,True]} in this case.

NDSolveValue::indexss: The DAE solver failed at t = 0.`. The solver is intended for index 1 DAE systems and structural analysis indicates that the DAE is structurally singular. Code here

Answer 1

We use scaling as in our model on this page.See also Vahl Davis, G.de (1983) : Natural convection of air in a square cavity : A bench mark numerical solution. Int. J. Numer. Methods Fluids 3, 249-264.
A system of equations describing free convection can be written in nondimensional form as follows $\nabla .\vec u=0, \frac{d\vec u}{dt}=Pr\nabla ^2\vec u -RaPrT\frac {\vec g}{g}$

$\frac {dT}{dt}=\nabla^2T$

$ d/dt =\partial/\partial t +(\vec u .\nabla )$

$\vec u$ is velocity field vector, $T$ - temperature, $\vec g $ - gravity vector, $Pr$ - the Prandtl number, $Ra$ - the Rayleigh number.

Note, that the temperature scale 40-(-40)=80 is included in Ra definition. Code should be modified as follows

Needs["NDSolve`FEM`"];
ResourceFunction["FEMAddOnsInstall"][];
Needs["FEMAddOns`"]
ClearAll[length, height, offsetx1, offsetx2, offsety];
sizes = {length -> 21/10, height -> 15/10,
   offsetx1 -> 1/10, offsetx2 -> 1/10,
   offsety1 -> 1/10, offsety2 -> 2/10
   };

Geometry
\[CapitalOmega]1 = Rectangle[{0, 0}, {length, height}] /. sizes;
\[CapitalOmega]2 = 
  Rectangle[{0, 0} + {0.2, 1.5 - 0.049}, {0.69, 0.049} + {0.2, 
     1.5 - 0.049}];

\[CapitalOmega]D1 = 
  RegionDifference[\[CapitalOmega]1, \[CapitalOmega]2];
\[CapitalOmega]3 = 
  Rectangle[{0, 0} + {1.2, 1.5 - 0.049}, {0.69, 0.049} + {1.2, 
     1.5 - 0.049}];
\[CapitalOmega]D2 = 
  RegionDifference[\[CapitalOmega]D1, \[CapitalOmega]3];
Mesh

\[CapitalOmega] = 
  DiscretizeRegion[\[CapitalOmega]D2, MaxCellMeasure -> 0.001, 
   AccuracyGoal -> 5, PrecisionGoal -> 5];
meshes = ToBoundaryMesh /@ {\[CapitalOmega]1, \[CapitalOmega]2, \
\[CapitalOmega]3};
bmesh = FEMUtils`BoundaryElementMeshJoin @@ meshes;
boxCoordinate = {1.05, 0.8};
PCMCoordinate1 = {0.5, 1.46};
PCMCoordinate2 = {1.5, 1.46};
markerColors = {Blue, Orange, Green};
markerCoordinate = {{boxCoordinate}, {PCMCoordinate1}, \
{PCMCoordinate2}};
markerspec = MapIndexed[{First@#1, First@#2} &, markerCoordinate];

mesh = ToElementMesh[bmesh, "RegionMarker" -> markerspec, 
   "MaxCellMeasure" -> 5*10^-3];
mesh["Wireframe"[]]
mesh["Wireframe"[
  "MeshElementStyle" -> 
   Map[Directive[FaceForm[#], EdgeForm[]] &, markerColors]]]

Show[{bmesh["Wireframe"], 
  Graphics[MapThread[{PointSize[0.02], #1, 
      Point /@ #2} &, {markerColors, markerCoordinate}]]}]
Material Properties
region = <|"air" -> 1, "pcm" -> 2|>;
vars = {T[t, x, y], t, {x, y}};
HeatTransferPDEComponent[vars, <|"MassDensity" -> \[Rho], 
  "SpecificHeatCapacity" -> Cp, "ThermalConductivity" -> k|>]
airpara = {Subscript[\[Rho], air] -> 
    QuantityMagnitude[ThermodynamicData["Air", "Density"]], 
        Subscript[Cp, air] -> 
    QuantityMagnitude[
     ThermodynamicData["Air", "IsobaricHeatCapacity"]],
        Subscript[k, air] -> 
    QuantityMagnitude[
     ThermodynamicData["Air", "ThermalConductivity"]]};
HS30Npara = {Subscript[Cp, pcmS] -> 2.1, Subscript[Cp, pcmL] -> 2.7, 
        Subscript[k, pcmS] -> 2.34, Subscript[k, pcmL] -> 0.6,
        Subscript[\[Rho], pcmS] -> 1430, latent -> 200, 
        tpmin -> -34, tpmax -> -26 };
PUFpara = {\[Rho] -> 50, cp -> 1500, k -> 0.025};

pars = <||>;
pars["MassDensity"] = 
  If[ElementMarker == region["air"], Subscript[\[Rho], air], 
   Subscript[\[Rho], pcmS]];
pars["SpecificHeatCapacity"] = 
  If[ElementMarker == region["air"], Subscript[k, air], Subscript[k, 
   pcmS] ];
pars["ThermalConductivity"] = 
  If[ElementMarker == region["air"], Subscript[k, air], 
   Subscript[k, pcmS]*IdentityMatrix[2]];
op = HeatTransferPDEComponent[vars, pars]

Initial and Boundary Condition
Tpcm = 0;
Tair = 0;
ic = {u[0, x, y] == 0,
   v[0, x, y] == 0,
   p[0, x, y] == 0,
   T[0, x, y] == 
    Piecewise[{{Tpcm, ElementMarker == region["pcm"]}, {Tair, 
       ElementMarker == region["air"]}}]};

Subscript[\[CapitalGamma], convectiveLeft] = 
  HeatTransferValue[
   x == 0 , {T[t, x, y], t, {x, y}}, <||>, <|
    "HeatTransferCoefficient" -> 0.7, "AmbientTemperature" -> 1|>];
Subscript[\[CapitalGamma], convectiveRight] = 
  HeatTransferValue[
   x == length /. sizes, {T[t, x, y], t, {x, y}}, <||>, <|
    "HeatTransferCoefficient" -> 0.7, "AmbientTemperature" -> 1|>];
Subscript[\[CapitalGamma], convectiveTop] = 
  HeatTransferValue[
   y == height /. sizes, {T[t, x, y], t, {x, y}}, <||>, <|
    "HeatTransferCoefficient" -> 0.7, "AmbientTemperature" -> 1|>];
Subscript[\[CapitalGamma], convectiveBottom] == 
  HeatTransferValue[
   y == 0 /. sizes, {T[t, x, y], t, {x, y}}, <||>, <|
    "HeatTransferCoefficient" -> 0.7, "AmbientTemperature" -> 1|>];
bcstemperatures = {
  DirichletCondition[
    T[t, x, y] == 
     0, {x, y} \[Element] \[CapitalOmega]2 || {x, 
       y} \[Element] \[CapitalOmega]3];
  }

ClearAll[\[Nu], \[Epsilon], Pr, Ra];
g = 9.8;
\[Alpha] = 0.0034;
\[Nu] = 1.48*10^-5;
\[Beta] = 2.17*10^-5;
L = (4*2.1*1.5)/(2*(2.1 + 1.5));
Temp = 40;
Subscript[T, \[Infinity]] = 0;
Ra -> ((g*\[Alpha])/(\[Nu]*\[Beta]))*(Temp - 
     Subscript[T, \[Infinity]])*L^3;
parameters = {Pr -> 0.7, Ra -> 10^6};

op = {
    
    
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y]
     + Inactive[Div][(-Pr 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}]
     + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y], 
        
    
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y]
     + Inactive[Div][(-Pr 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}]
     + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y] - Ra Pr T[t, x, y],
        
    
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y] + 
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y],
        
    
\!\(\*SuperscriptBox[\(T\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y]
     + Inactive[Div][(-Inactive[Grad][T[t, x, y], {x, y}]), {x, y}]
     + {u[t, x, y], v[t, x, y]} . Inactive[Grad][T[t, x, y], {x, y}]
    
    } /. parameters;
walls = {
   DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, {x, y}]
   };
reference = DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0];
bcsflow = {walls, reference };

Solution over mesh with $Ra=10^6$ takes about 120 s

tmax = 1;
Monitor[
 AbsoluteTiming[
  {xVel, yVel, pressure, temperature} = NDSolveValue[
     Flatten@{op == {0, 0, 0, 
         Subscript[\[CapitalGamma], convectiveLeft] + 
          Subscript[\[CapitalGamma], convectiveRight] + 
          Subscript[\[CapitalGamma], convectiveTop]}, bcsflow, 
       bcstemperatures, ic}, {u, v, p, T}, {x, y} \[Element] mesh, {t,
       0, tmax},
     Method -> {"TimeIntegration" -> {"IDA", 
         "MaxDifferenceOrder" -> 2},
       "PDEDiscretization" -> {
         "MethodOfLines", "SpatialDiscretization" -> {
           "FiniteElement",
           "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}}
         }
       },
     EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])
     ];
  ],
 currentTime
 ]

Visualization

{StreamDensityPlot[
  {xVel[tmax, x, y], yVel[tmax, x, y]}, {x, y} \[Element] 
   xVel["ElementMesh"],
  AspectRatio -> Automatic, PlotRange -> All,
  StreamPoints -> Fine, PlotLegends -> Automatic
  ], DensityPlot[
  temperature[tmax, x, y], {x, y} \[Element] 
   temperature["ElementMesh"],
  AspectRatio -> Automatic, PlotRange -> All,
  ColorFunction -> "TemperatureMap", PlotLegends -> Automatic
  ]}

Solution on $\Omega$ with $Ra=10^5$ takes about 90 s