# NDSolveValue taking too long

I am trying to solve heat diffusion equation on hollow cylinder with constant DirichletCondition on inner radius and zero NeumannCondition (no conduction) on outer radius. Initial condition is linearly with r decreasing temperature. My code:

h = 10;
cyl1 = Cylinder[{{0, 0, 0}, {0, 0, h}}, 2];
cyl2 = Cylinder[{{0, 0, 0}, {0, 0, h}}, 4];
cyl = DiscretizeRegion[RegionDifference[cyl2, cyl1]];
Laplacian[u[t, r, θ, z], {r, θ, z}, "Cylindrical"]
ClearAll[r]
ifun = NDSolveValue[{D[u[t, r, θ, z], t] ==
Laplacian[u[t, r, θ, z], {r, θ, z},
"Cylindrical"] + NeumannValue[0.1, (0 < z < h) && (r > 3.9)],
DirichletCondition[u[t, r, θ, z] == 50, r < 2.1],
u[0.001, r, θ, z] == 50 - ((50 - 10)/(4 - 2))*r},
u, {t, 1, 5}, {r, θ, z} ∈ cyl, Method -> Automatic];
SliceContourPlot3D[
ifun[4, r, θ, z], {r, θ, z} ∈ cyl,
ColorFunction -> "TemperatureMap", Boxed -> False, Axes -> None]


The question is why is NDSolveValue taking so long.

EDIT: How to convince Mathematica to interpret (r,theta,z) as cylindrical (not cartesian) coordinates of a region cyl?

• Hi kular. On your machine, how long does it take to run? Commented Apr 24, 2020 at 20:02
• I'm running it for 2 hours, still no result. Commented Apr 24, 2020 at 20:07
• @kular You made a typo with Cylinder[]. In the cylindrical coordinate it is not cylinder, but region $2\le r\le 2, 0\le \theta \le 2\pi, 0\le z\le h$ Commented Apr 25, 2020 at 19:21

Try

h = 10;
innerR = 2;
outerR = 4;
cyl1 = Cylinder[{{0, 0, 0}, {0, 0, h}}, innerR];
cyl2 = Cylinder[{{0, 0, 0}, {0, 0, h}}, outerR];
cyl = ToElementMesh[RegionDifference[cyl2, cyl1], MaxCellMeasure -> .5];
cyl["Wireframe"]


pde = D[u[t, r, θ, z], t] == Laplacian[u[t, r, θ, z], {r, θ, z}, "Cylindrical"] +
NeumannValue[1/10, (0 < z < h) && (r > 39/10)];

ic = u[1/1000, r, θ, z] == 50 - ((50 - 10)/(4 - 2))*r;

ifun = NDSolveValue[{pde, DirichletCondition[u[t, r, θ, z] == 50,
r < 21/10], ic}, u, {t, 1/1000, 5}, {r, θ, z} ∈ cyl]


Finishes right away.

• There is a typo with Cylinder[]. In the cylindrical coordinates it is not cylinder, but Coboud[] region - see Tim Laska answer and my also. Commented Apr 25, 2020 at 23:02

In cylindrical space, the annular region would look like a cuboid. Since there is no $$\theta$$ dependence in your model, you may as well do an axisymmetric model. Fortunately, Mathematica has a heat transfer verification manual to give a guide on how to setup a verified case, from which, I will shamelessly copy and paste.

ClearAll[HeatTransferModel]
HeatTransferModel[T_, X_List, k_, \[Rho]_, Cp_, Velocity_, Source_] :=
Module[{V, Q, a = k},
V = If[Velocity === "NoFlow",
Q = If[Source === "NoSource", 0, Source];
If[ FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
If[ VectorQ[a], a = DiagonalMatrix[a]];
a = PiecewiseExpand[Piecewise[{{-a, True}}]];
Inactive[Div][a.Inactive[Grad][T, X], X] + V - Q]

HeatTransferModelAxisymmetric[T_, {r_, z_}, k_, \[Rho]_, Cp_,
Velocity_, Source_] :=
Module[{V, Q},
V = If[Velocity === "NoFlow",
Q = If[Source === "NoSource", 0, Source];
1/r*D[-k*r*D[T, r], r] + D[-k*D[T, z], z] + V - Q]

TimeHeatTransferModel[T_, TimeVar_, X_List, k_, \[Rho]_, Cp_,
Velocity_, Source_] := \[Rho]*Cp*D[T, {TimeVar, 1}] +
HeatTransferModel[T, X, k, \[Rho], Cp, Velocity, Source]

TimeHeatTransferModelAxisymmetric[T_, TimeVar_, {r_, z_}, k_, \[Rho]_,
Cp_, Velocity_, Source_] :=
\[Rho]*Cp*D[T, {TimeVar, 1}] +
HeatTransferModelAxisymmetric[T, {r, z}, k, \[Rho], Cp, Velocity,
Source]


Now, we can setup a verified heat equation operator for an axisymmetric case and solve your problem.

h = 10;
rmin = 2;
rmax = 4;
Subscript[\[CapitalGamma], temp] =
DirichletCondition[u[t, r, z] == 50, r == rmin];
nv = NeumannValue[0.1, r == rmax];
ic = {u[0, r, z] == 50 - ((50 - 10)/(4 - 2))*r};
\[CapitalOmega] = Rectangle[{rmin, 0}, {rmax, h}];
tend = 5;
parmop = TimeHeatTransferModelAxisymmetric[u[t, r, z], t, {r, z},
k, \[Rho], Cp, "NoFlow", "NoSource"];
op = parmop /. {k -> 1, \[Rho] -> 1, Cp -> 1};
pde = {op == nv, Subscript[\[CapitalGamma], temp], ic};
ifun = NDSolveValue[pde,
u, {t, 0, tend}, {r, z} \[Element] \[CapitalOmega]]


The solution returns quickly. We can visualize the solution using the example in the verification manual or from the example on the Wolfram website shown here.

uRange = MinMax[ifun["ValuesOnGrid"]];
legendBar =
BarLegend[{"TemperatureMap", uRange}, 50,
LegendLabel -> Style["[\[Degree]C]", Opacity[0.6]]];
options = {PlotRange -> uRange,
ColorFunction -> ColorData[{"TemperatureMap", uRange}],
ContourStyle -> Opacity[0.1], ColorFunctionScaling -> False,
Contours -> 30, AspectRatio -> 1,
PlotPoints -> 41, FrameLabel -> {"r", "z"},
PlotLabel -> Style["Temperature Field: u(t,r,z)", 18],
AspectRatio -> Automatic, ImageSize -> 250};
nframes = 80;
frames = Table[
Legended[
ContourPlot[ifun[t, r, z], {r, z} \[Element] \[CapitalOmega],
Evaluate[options]], legendBar], {t, 0, tend, tend/nframes}];
frames = (Rasterize[#1, "Image", ImageResolution -> 80] &) /@ frames;
ListAnimate[frames, SaveDefinitions -> True]
cpfn = ContourPlot[ifun[#, r, z], {r, z} \[Element] \[CapitalOmega],
ColorFunction -> "Temperature"] &;
frames2 =
Rasterize[#1, "Image",
ImageResolution ->
100] & /@ (Show[
RegionPlot3D[
rmin^2 <= x^2 + y^2 <= rmax^2 && 0 <= z <= h, {x, -rmax,
rmax}, {y, -rmax, rmax}, {z, 0, h}, Boxed -> False,
Axes -> False, PlotPoints -> 40, PlotStyle -> {Opacity[0.2]},
Mesh -> False],
Graphics3D[{EdgeForm[Red], FaceForm[Gray],
GraphicsComplex[{{rmin, 0, 0}, {rmax, 0, 0}, {rmax, 0,
h}, {rmin, 0, h}}, {Texture[
Show[cpfn[#], Frame -> False, PlotRangePadding -> None]],
Lighting -> {{"Ambient", White}},
Polygon[{{1, 2, 3, 4}},
VertexTextureCoordinates -> {{{0, 0}, {1, 0}, {1, 1}, {0,
1}}}]}]}], ImageSize -> 200] & /@
Subdivide[0, tend, 80]);
ListAnimate[frames2, SaveDefinitions -> True]


Here is an example of how to make it look more 3D with a SliceContourPlot3D

cyl1 = Cylinder[{{0, 0, 0}, {0, 0, h}}, 2];
cyl2 = Cylinder[{{0, 0, 0}, {0, 0, h}}, 4];
cyl = RegionDifference[cyl2, cyl1];
frames3 =
Rasterize@
SliceContourPlot3D[
ifun[#, Sqrt[x^2 + y^2], z], {x, y, z} \[Element] cyl,
PlotRange -> uRange,
ColorFunction -> ColorData[{"TemperatureMap", uRange}],
ContourStyle -> Opacity[0.5], ColorFunctionScaling -> False,
Contours -> 30, Boxed -> False, Axes -> False,
PlotPoints -> 40] & /@ Subdivide[0, tend, 40];
ListAnimate@frames3


• (+1) But this is actually 2D, not 3D. How it could be solve in a case of 3D? Commented Apr 25, 2020 at 19:07
• @AlexTrounev Thanks! Yes, it is an axisymmetric 2D model. The OP question didn't have any $\theta$ dependence, so why not take advantage of symmetry? I was having difficulty coming up with an advantage to the approach in the OP unless it is to avoid some simple post processing steps. Your answer revealed something subtle, but important. Increasing the number of boundary surfaces going from cartesian (4) to cylindrical (6), increases the risk of boundary-boundary artefacts. Those can be difficult to fix. You may need to control PointElement markers to fix. Commented Apr 25, 2020 at 19:47
• I am wonder that it working in 3D with periodic boundary condition. Usually I solve this problem in cartesian coordinates, then I have no headache with boundary elements and PeriodicBoundaryCondition[]. Commented Apr 25, 2020 at 22:54

In cylindrical coordinates region Cylinder[] should be transformed into Cuboid[], so effectively we should solve this problem in the region {r, 2, 4}, {\[Theta], 0, 2 Pi}, {z, 0, h} with a periodic boundary condition on $$\theta$$. Code

Needs["NDSolveFEM"]; h = 10; reg =
ImplicitRegion[
2 <= r <= 4 && 0 <= \[Theta] <= 2 Pi && 0 <= z <= h, {r, \[Theta],
z}];

mesh = ToElementMesh[reg]

mesh["Wireframe"]
pbc = PeriodicBoundaryCondition[u[t, r, \[Theta], z], \[Theta] == 0,
TranslationTransform[{0, 2 Pi, 0}]];

ifun = NDSolveValue[{D[u[t, r, \[Theta], z], t] -
Laplacian[u[t, r, \[Theta], z], {r, \[Theta], z},
"Cylindrical"] == NeumannValue[0., True],
DirichletCondition[u[t, r, \[Theta], z] == 50,
r < 2.1 && 0 < \[Theta] < 2 Pi],
u[1, r, \[Theta], z] == 50 - ((50 - 10)/(4 - 2))*r},
u, {t, 1, 5}, {r, \[Theta], z} \[Element] mesh]

SliceContourPlot3D[
ifun[4, r, \[Theta],
z], "ZStackedPlanes", {r, \[Theta], z} \[Element] mesh,
ColorFunction -> "TemperatureMap", Boxed -> False,
AxesLabel -> Automatic]


Projection into {x,y,z} coordinates:

SliceContourPlot3D[
ifun[4, Sqrt[x^2 + y^2], ArcTan[x, y] + Pi,
z], "ZStackedPlanes", {x, -4, 4}, {y, -4, 4}, {z, 0, h},
ColorFunction -> "TemperatureMap", Boxed -> False,
AxesLabel -> Automatic]


• Should there be any $\theta$ dependence? I would have thought there would be a series of straight lines. Any thoughts? Commented Apr 25, 2020 at 17:56
• @TimLaska Probably we should increase number of points on $\theta$ or may be it just issues of PeriodicBoundaryCondition. But I see that you post a right model. Commented Apr 25, 2020 at 18:41
• I think there is an issue where you have two types nodal boundary conditions on adjacent edges competing for the same corner point. Mathematica makes a hidden decision about which one takes precedence. Resolution will help mitigate the issue, but it would be nice to change the precedence so that you could experiment. Commented Apr 25, 2020 at 18:56
• @TimLaska It was issues of Method->Automatic`. With FEM it looks better. Commented Apr 25, 2020 at 20:10
• Indeed, that looks much better! Commented Apr 25, 2020 at 20:30