NDSolveValue taking too long

Question

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?

Answer 1

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.

Answer 2

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", 
    0, \[Rho]*Cp*Velocity.Inactive[Grad][T, X]];
  Q = If[Source === "NoSource", 0, Source];
  If[ FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
  If[ VectorQ[a], a = DiagonalMatrix[a]];
  (* Note the - sign in the operator *)
  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", 
    0, \[Rho]*Cp*Velocity.Inactive[Grad][T, {r, z}]];
  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

Answer 3

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["NDSolve`FEM`"]; 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]