# Speeding up Boussinesq equations solving?

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

• Crossposted here. Commented May 20, 2022 at 21:24
• @RohitNamjoshi yes, it's the same person. Do you have any suggestions? Commented May 20, 2022 at 21:29
• I do not. When you crosspost the same question on multiple sites it is good etiquette to add crosslinks. That way the people trying to help you can check if the question has already been answered elsewhere and if so they do not have to waste time answering it. Commented May 20, 2022 at 23:16
• @LionSahara Please, check, that the Rayleigh number in your problem is about 2.22412*10^10, and therefore this is turbulent convection. Also mesh you used is very rough for this problem. Even for laminar convection with $Ra=10^5$ I have used option "MaxCellMeasure" -> 0.001 to compare with other solvers - see my post on community.wolfram.com/groups/-/m/t/1433064?p_p_auth=KT89yn6t Commented May 21, 2022 at 4:35
• @AlexTrounev Thank you for your suggestion. I changed Ra value and Maxcellmeasure as you suggested. The simulation time took forever. The simulation has been running for an hour and t is roughly 0.26 steps. I have no idea what is the problem here. Commented May 21, 2022 at 12:10

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["NDSolveFEM"];
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 = FEMUtilsBoundaryElementMeshJoin @@ 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"],
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

• Thank you very much. I get my code works perfectly. Thank you again. Commented May 22, 2022 at 13:28
• You are welcome! Commented May 22, 2022 at 15:11
• I have a question. I understand that the temperature scale 40-(-40)=80. Why do we have to set up the initial temperature of PCM at 0 and Tatm is 1? What if I want to use Tpcm at -25 and Tatm at 30. Do I need to calculate a new Ra? Commented May 22, 2022 at 19:57
• Yes, any temperature scale should be included in Ra definition, since force depends on temperature scale included in Ra. Commented May 22, 2022 at 23:51
• Thank you very much. I understand the equation like a crystal. Commented May 23, 2022 at 11:43