# How to diminish computation time when nonlinearity appears in 2D heat conduction equation?

I am trying to simulate heating and melting of the steel plate by means of FEM.The model is based on nonlinear heat conduction equation in axial symmetry case.

The problem statement is the next: $$\rho c_{eff}\frac{\partial T}{\partial t}= \frac{1}{r}\frac{\partial}{\partial r}\left(r\lambda \frac{\partial T}{\partial r} \right) + \frac{\partial}{\partial z}\left(\lambda \frac{\partial T}{\partial z} \right),\\ 0\leq r\leq L_r,~0\leq z\leq L_z,~0\leq t\leq t_f$$ $$\lambda \frac{\partial T}{\partial z}\Bigg|_{z=L_z}=q_{0}exp(-a r^2),~~\frac{\partial T}{\partial r}\Bigg|_{r=L_r}=0, T|_{z=0}=T_0\\T(0,r,z)=T_0$$

To take into account latent heat of fusion $$L$$ the effective heat capacity is introduced $$c_{eff}=c_{s}(1-\phi)+c_{l}\phi+ L\frac{d \phi}{dT}$$, where $$\phi$$ is a fraction of liquid phase, $$c_s, c_l$$ are the heat capacity of solid and liquid phase respectively. Smoothed Heaviside function

$$h(x,\delta)=\left\{\begin{array}{l l l} 0,& x<-\delta\\ 0.5\left(1+\frac{x}{\delta}+\frac{1}{\pi}sin(\frac{\pi x}{\delta}) \right), &\mid x\mid\leq \delta\\ 1,& x>\delta \end{array} \right.$$

is employed to describe mushy zone so that $$\phi(T)=h(T-T_m,\Delta T_{m}/2)$$, where $$T_m$$ and $$\Delta T_m$$ are melting temperature and melting range respectively. FE approximation is used for spatial discretization of PDE whereas time derivative is approximated by means of first order finite difference scheme: $$\left.\frac{\partial T}{\partial t}\right|_{t=t^{k}} \approx \frac{T(t^k,r,z)-T(t^{k-1},r,z)}{\tau}$$

where $$\tau$$ is a time step size. For calculation of $$c_{eff}$$ at k-th time step the temperature field from k-1 time step is utilized. After discretization in time one can rewrite equation:

$$c_{eff}\left(T(t^{k-1},r,z)\right) \frac{T(t^k,r,z)-T(t^{k-1},r,z)}{\tau}=\frac{1}{r}\frac{\partial}{\partial r}\left(r\lambda \frac{\partial T(t^k,r,z)}{\partial r} \right) + \frac{\partial}{\partial z}\left(\lambda \frac{\partial T(t^k,r,z)}{\partial z} \right)$$

At each time step the DampingCoefficients is corrected in InitializePDECoefficients[] so that interpolation is used for $$c_{eff}$$.Such approach leads to significant grows of computational time in comparison with solution of linear problem when $$c_{eff}$$=const. I also tried to use ElementMarker to set certain value of $$c_{eff}$$ in each element. Such approach allows to avoid interpolation but the computation time is getting more larger. This last fact I can not understand at all. As to me the duration of FE matrix assembly should be diminished when interpolation for $$c_{eff}$$ is avoided.

Needs["NDSolveFEM"];
Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];


Setting of the computational domain dimensions and mesh generation:

Lr = 2*10^-2; (*dimension of computational domain in r-direction*)
Lz = 10^-2;   (*dimension of computational domain in z-direction*)
mesh = ToElementMesh[FullRegion, {{0, Lr}, {0, Lz}},MaxCellMeasure -> {"Length" -> Lr/50}, "MeshOrder" -> 1]
mesh["Wireframe"] Input parameters of the model:

lambda = 22;         (*heat conductivity*)
density = 7200;      (*density*)
Cs = 700;            (*specific heat capacity of solid*)
Cl = 780;            (*specific heat capacity of liquid*)
LatHeat = 272*10^3;  (*latent heat of fusion*)
Tliq = 1812;         (*melting temperature*)
MeltRange = 100;     (*melting range*)
To = 300;            (*initial temperature*)
SPow = 1000;         (*source power*)
R = Lr/4;            (*radius of heat source spot*)
a = Log/R^2;
qo = (SPow*a)/Pi;
q[r_] := qo*Exp[-r^2*a]; (*heat flux distribution*)
tau = 10^-3;         (*time step size*)
ProcDur = 0.2;       (*process duration*)


Smoothed Heaviside function:

Heviside[x_, delta_] := Module[{res},
res = Piecewise[

{
{0, Abs[x] < -delta},
{0.5*(1 + x/delta +  1/Pi*Sin[(Pi*x)/delta]), Abs[x] <= delta},
{1, x > delta}
}

];
res
]


Smoothed Heaviside function derivative:

HevisideDeriv[x_, delta_] := Module[{res},
res = Piecewise[
{

{0, Abs[x] > delta},

{1/(2*delta)*(1 + Cos[(Pi*x)/delta]), Abs[x] <= delta}
}
];
res
]


Effective heat capacity:

EffectHeatCapac[tempr_] := Module[{phase},
phase = Heviside[tempr - Tliq, MeltRange/2];
Cs*(1 - phase) + Cl*phase +LatHeat*HevisideDeriv[tempr - Tliq, 0.5*MeltRange]
]


Numerical solution of PDE:

ts = AbsoluteTime[];

vd = NDSolveVariableData[{"DependentVariables" -> {u},"Space" -> {r,z},"Time" -> t}];
sd = NDSolveSolutionData[{"Space","Time"} -> {ToNumericalRegion[mesh], 0.}];

DirichCond=DirichletCondition[u[t, r, z] ==To,z==0];
NeumCond=NeumannValue[q[r],z==Lz];
initBCs=InitializeBoundaryConditions[vd,sd, {{DirichCond, NeumCond}}];
methodData = InitializePDEMethodData[vd, sd] ;
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

xlast = Table[{To}, {methodData["DegreesOfFreedom"]}];
TemprField = ElementMeshInterpolation[{mesh}, xlast];
NumTimeStep = Floor[ProcDur/tau];

For[i = 1, i <= NumTimeStep, i++,

(*
(*Setting of PDE coefficients for linear problem*)
pdeCoefficients=InitializePDECoefficients[vd,sd,"ConvectionCoefficients"->     {{{{-lambda/r, 0}}}},
"DiffusionCoefficients" -> {{-lambda*IdentityMatrix}},
"DampingCoefficients" -> {{Cs*density}}];
*)

(*Setting of PDE coefficients for nonlinear problem*)

pdeCoefficients =
InitializePDECoefficients[vd, sd,
"ConvectionCoefficients" -> {{   {{-(lambda/r), 0}}  }},
"DiffusionCoefficients" -> {{-lambda*IdentityMatrix}},
"DampingCoefficients" -> {{EffectHeatCapac[TemprField[r, z]]*
density}}];

discretePDE = DiscretizePDE[pdeCoefficients, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];
discreteBCs];

A = damping/tau + stiffness;

x = LinearSolve[A,b,Method -> {"Krylov", Method -> "BiCGSTAB",
"Preconditioner" -> "ILU0","StartingVector"->Flatten[xlast,1]}];
TemprField = ElementMeshInterpolation[{mesh}, x];
xlast = x;
]
te = AbsoluteTime[];
te - ts


Visualization of the calculation results

ContourPlot[TemprField[r, z], {r, z} \[Element] mesh,
AspectRatio -> Lz/Lr, ColorFunction -> "TemperatureMap",
Contours -> 50, PlotRange -> All,
PlotLegends -> Placed[Automatic, After], FrameLabel -> {"r", "z"},
PlotPoints -> 50, PlotLabel -> "Temperature field", BaseStyle -> 16] On my laptop the computation time are 63 sec and 2.17 sec for nonlinear and linear problems respectively.This question can be generalized to the case when $$\lambda=\lambda(T)$$. I would appreciate if anyone could please show me a good way which leads to time savings. Thanks in advance for your help.

• The code seems to be incomplete. For instance, methodData is undefined. Please edit your post and add all relevant code. Mar 21, 2019 at 9:23
• When you use the method "Krylov" for transient PDE, it is usually a very good idea to supply the solution of the linear system of the previous time iteration as value for "StartingVector". You should also try to find suitable values for the "Tolerance" option; the default values are usually way too low. Moreover, using the backend SparseArrayKrylovLinearSolve directly is usually a bit faster and also allows for reusing the preconditioner (generated by SparseArraySparseMatrixILU and applied with SparseArraySparseMatrixApplyILU). Mar 21, 2019 at 9:30
• When you will have added the rest of the code I may show you how to use SparseArrayKrylovLinearSolve; so please ping me with @Henrik in a comment when you are done. Mar 21, 2019 at 9:32
• Well, then a complete code would help. I have to test against something... Mar 21, 2019 at 12:26
• I would be good if you could give a complete code, without it is just guessing. You could also split the InitilizePDECoefficient into a linear and a nonlinear part. Discretize the linear part before the loop and then do linearDiscretePDE["SystemMatrices"] + nonlinearDiscretePDE["SystemMatrices"]. Before using "Krylov" try Method->"Pardiso" for LinearSolve. Mar 22, 2019 at 5:53

As promised, here my 6 pence.

### Basic settings

Needs["NDSolveFEM"];
Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];

Lr = 2*10^-2;(*dimension of computational domain in r-direction*)
Lz = 10^-2;(*dimension of computational domain in z-direction*)
mesh = ToElementMesh[FullRegion, {{0, Lr}, {0, Lz}}, MaxCellMeasure -> {"Length" -> Lr/50}, "MeshOrder" -> 1]
mesh["Wireframe"]

lambda = 22.;         (*heat conductivity*)
density = 7200.;      (*density*)
Cs = 700.;            (*specific heat capacity of solid*)
Cl = 780.;            (*specific heat capacity of liquid*)
LatHeat = 272.*10^3;  (*latent heat of fusion*)
Tliq = 1812.;         (*melting temperature*)
MeltRange = 100.;     (*melting range*)
To = 300.;            (*initial temperature*)
SPow = 1000.;         (*source power*)
R = Lr/4.;            (*radius of heat source spot*)
a = Log[100.]/R^2;
qo = (SPow*a)/Pi;
q[r_] := qo*Exp[-r^2*a]; (*heat flux distribution*)
tau = 10^-3;         (*time step size*)
ProcDur = 0.2;       (*process duration*)

Heviside[x_, delta_] := Piecewise[{{0,
Abs[x] < -delta}, {0.5*(1 + x/delta + 1/Pi*Sin[(Pi*x)/delta]),
Abs[x] <= delta}, {1, x > delta}}];

HevisideDeriv[x_, delta_] := Piecewise[{{0,
Abs[x] > delta}, {1/(2*delta)*(1 + Cos[(Pi*x)/delta]),
Abs[x] <= delta}}];

EffectHeatCapac[tempr_] := Module[{phase},
phase = Heviside[tempr - Tliq, MeltRange/2];
Cs*(1 - phase) + Cl*phase + LatHeat*HevisideDeriv[tempr - Tliq, 0.5*MeltRange]];


### Compiled versions of smoothed Heaviside functions

cHeaviside = Compile[{{x, _Real}, {delta, _Real}},
Piecewise[{
{0.,
Abs[x] < -delta}, {0.5*(1 + x/delta + 1./Pi*Sin[(Pi*x)/delta]),
Abs[x] <= delta}, {1., x > delta}}
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
];
cHeavisideDeriv = Compile[{{x, _Real}, {delta, _Real}},
Piecewise[{
{0., Abs[x] > delta},
{1./(2*delta)*(1. + Cos[(Pi*x)/delta]), Abs[x] <= delta}}
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
];
cEffectHeatCapac[tempr_] :=
With[{phase = cHeaviside[tempr - Tliq, MeltRange/2]},
Cs*(1 - phase) + Cl*phase + LatHeat*cHeavisideDeriv[tempr - Tliq, 0.5*MeltRange]
];


### A Fast matrix asssembler routine

Copied from here.

SetAttributes[AssemblyFunction, HoldAll];

Assembly::expected = "Values list has 2 elements. Expected are 1 elements. Returning  prototype.";

Assemble[pat_?MatrixQ, dims_, background_: 0.] :=
Module[{pa, c, ci, rp, pos},
pa = SparseArraySparseArraySort@SparseArray[pat -> _, dims];
rp = pa["RowPointers"];
ci = pa["ColumnIndices"];
c = Length[ci];
pos = cLookupAssemblyPositions[Range[c], rp, Flatten[ci], pat];
Module[{a},
a = <|
"Dimensions" -> dims,
"Positions" -> pos,
"RowPointers" -> rp,
"ColumnIndices" -> ci,
"Background" -> background,
"Length" -> c
|>;
AssemblyFunction @@ {a}]
];

AssemblyFunction /: a_AssemblyFunction[vals0_] :=
Module[{len, expected, dims, u, vals, dat},
dat = a[];
If[VectorQ[vals0], vals = vals0, vals = Flatten[vals0]];
len = Length[vals];
expected = Length[dat[["Positions"]]];
dims = dat[["Dimensions"]];
If[len === expected,
If[Length[dims] == 1,
u = ConstantArray[0., dims[]];
u[[dat[["ColumnIndices"]]]] = AssembleDenseVector[dat[["Positions"]], vals, {dat[["Length"]]}];
u,
SparseArray @@ {Automatic, dims,
dat[["Background"]], {1, {dat[["RowPointers"]],
dat[["ColumnIndices"]]},
AssembleDenseVector[dat[["Positions"]],
vals, {dat[["Length"]]}]}}
],
Message[Assembly::expected, expected, len];
Abort[]]
];

cLookupAssemblyPositions =
Compile[{{vals, _Integer, 1}, {rp, _Integer, 1}, {ci, _Integer, 1}, {pat, _Integer, 1}},
Block[{k, c, i, j},
i = CompileGetElement[pat, 1];
j = CompileGetElement[pat, 2];
k = CompileGetElement[rp, i] + 1;
c = CompileGetElement[rp, i + 1];
While[k < c + 1 && CompileGetElement[ci, k] != j,
++k
];
CompileGetElement[vals, k]
],
RuntimeAttributes -> {Listable},
Parallelization -> True,
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];

AssembleDenseVector =
Compile[{{ilist, _Integer, 1}, {values, _Real, 1}, {dims, _Integer, 1}}, Block[{A}, A = Table[0., {CompileGetElement[dims, 1]}];
Do[A[[CompileGetElement[ilist, i]]] +=
CompileGetElement[values, i], {i, 1, Length[values]}];
A
],
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];


### Damping matrix assembly code

Mostly reverse engineered, so I am actually not 100% sure that this does what it should...

As far a I got it, the damping matrix with respect to function $$f \colon \varOmega \to \mathbb{R}$$ should encode the bilinear form

$$(u,v) \mapsto \int_{\varOmega} u(x) \, v(x) \, f(x) \, \mathrm{d} x.$$ in terms of the FEM basis functions. Since the FEM basis functions have very local support, we go over the finite elements of the mesh (quads in this case) and compute the local contributions to the overall bilinear form. Then it is a matter of index juggling to assemble the

This assumes bi-linear interpolation on quads and employs Gaussian quadrature with 2 integration points per dimension for integration. (For triangular or tetrahedral meshes, exact integration can be used instead.)

(* for each quad, getWeakLaplaceCombinatoricsQuad is supposed to produce the $$i-j$$-indices of each of the 16 entries of the local $$4 \times 4$$ metrix within the global matrix.*)
With[{code = Flatten[Table[Table[{
CompileGetElement[q, i],
CompileGetElement[q, j]
}, {i, 1, 4}], {j, 1, 4}], 1]},
Compile[{{q, _Integer, 1}},
code,
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]
]
];

(* this snippet computes the symbolic expression for the local matrices and then compiles it into the function getLocalDampingMatrices*)
Block[{dim, PP, UU, FF, p, u, f, integrant, x, ω, localmatrix},
dim = 2;
PP = Table[CompileGetElement[P, i, j], {i, 1, 4}, {j, 1, dim}];
UU = Table[CompileGetElement[U, i], {i, 1, 4}];
FF = Table[CompileGetElement[F, i], {i, 1, 4}];

(* bi-linear interpolation of the quadrilateral; maps the standard quare onto the quadrilateral defined by PP[], PP[], PP[], PP[]*)
p = {s, t} \[Function] (PP[] (1 - s) + s PP[]) (1 - t) + t (PP[] (1 - s) + s PP[]);
(* bi-linear interpolation of the function values of u*)
u = {s, t} \[Function] (UU[] (1 - s) + s UU[]) (1 - t) + t (UU[] (1 - s) + s UU[]);
(* bi-linear interpolation of the function values of f*)
f = {s, t} \[Function] (FF[] (1 - s) + s FF[]) (1 - t) + t integrant = {s, t} \[Function] Evaluate[f[s, t] u[s, t]^2 Abs[Det[D[p[s, t], {{s, t}, 1}]]]];
{x, ω} = Most[NIntegrateGaussRuleData[2, MachinePrecision]];

(* using D to obtain the local matrix from its quadratic form*)
localmatrix = 1/2 D[
Flatten[KroneckerProduct[ω, ω]].integrant @@@ Tuples[x, 2],
{UU, 2}
];

(* getLocalDampingMatrices computes the local $$4 \times 4$$-matrices from the quad vertex coordinates P (supposed to be a $$4 \times 2$$-matrix) and from the function values F (supposed to be a $$4$$-vector) *)
getLocalDampingMatrices = With[{code = localmatrix},
Compile[{{P, _Real, 2}, {F, _Real, 1}},
code,
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]
];
];

Module[{fdata, localmatrices},
assembler[Flatten[localmatrices]]
];


The function getDampingMatrix eats an AssemblyFunction object assembler_, the list quads of of all quads (as a list of 4-vectors of the vertex indices), the list quaddata (a list of $$4 \times 2$$-matrix with the vertex positions, and a list fvals with the values of the function $$f$$ at the vertices of the mesh. It spits out the completely assembled damping matrix.

### Using DiscretizePDE only once

This requires the old implementation of EffectHeatCapac.

u =.
vd = NDSolveVariableData[{"DependentVariables" -> {u}, "Space" -> {r, z}, "Time" -> t}];
sd = NDSolveSolutionData[{"Space", "Time"} -> {ToNumericalRegion[mesh], 0.}];

DirichCond = DirichletCondition[u[t, r, z] == To, z == 0];
NeumCond = NeumannValue[q[r], z == Lz];
initBCs = InitializeBoundaryConditions[vd, sd, {{DirichCond, NeumCond}}];
methodData = InitializePDEMethodData[vd, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

x0 = ConstantArray[To, {methodData["DegreesOfFreedom"]}];
TemprField = ElementMeshInterpolation[{mesh}, x0];
NumTimeStep = Floor[ProcDur/tau];

pdeCoefficients = InitializePDECoefficients[vd, sd,
"ConvectionCoefficients" -> {{{{-(lambda/r), 0}}}},
"DiffusionCoefficients" -> {{-lambda*IdentityMatrix}},
"DampingCoefficients" -> {{EffectHeatCapac[TemprField[r, z]] density}}
];
discretePDE = DiscretizePDE[pdeCoefficients, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];


### Running the simulation

By removing the bottlenecks DiscretizePDE and (much more severely) ElementMeshInterpolation, the loop requires now only 0.32 seconds to execute. We also profit from the fact that, by utilizing the AssemblyFunction assembler, we don't have to recompute any sparse array patterns . Moreover, utilizing an undocumented syntax for the SparseArray constructor circumvents certain further performance degradations .

So this is now faster by a factor of 100.

x = x0;
taustiffness = tau stiffness;

Do[
A = damping + taustiffness;
x = LinearSolve[A, b, Method -> {"Krylov",
Method -> "BiCGSTAB",
"Preconditioner" -> "ILU0",
"StartingVector" -> x
}
];
,
{i, 1, NumTimeStep}]; // AbsoluteTiming // First


0.325719

### Using ElementMeshInterpolation only once on the end for plotting

TemprField = ElementMeshInterpolation[{mesh}, x];

ContourPlot[TemprField[r, z], {r, z} ∈ mesh,
AspectRatio -> Lz/Lr,
ColorFunction -> "TemperatureMap",
Contours -> 50,
PlotRange -> All,
PlotLegends -> Placed[Automatic, After],
FrameLabel -> {"r", "z"},
PlotPoints -> 50,
PlotLabel -> "Temperature field",
BaseStyle -> 16] After running

fvals = cEffectHeatCapac[x] density;


the line

assembler[localmatrices];


is basically equivalent to using SparseArray for additive assembly as follows:

(* switching to additive matrix assembly *)
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}];
SparseArray[pat -> Flatten[localmatrices], {1, 1} Length[fvals], 0.];


Maybe this helps to understand how getWeakLaplaceCombinatoricsQuad and getLocalDampingMatrices are related.

I implemented a somewhat slicker interface for simplicial meshes of arbitrary dimensions here.

So let's assume that we started with the following triangle mesh:

mesh = ToElementMesh[FullRegion, {{0, Lr}, {0, Lz}},
MaxCellMeasure -> {"Length" -> Lr/50}, "MeshOrder" -> 1,
MeshElementType -> TriangleElement];


Then one has to convert the mesh once into a MeshRegion.

Ω = MeshRegion[mesh];


damping = getDampingMatrix[assembler, quads, quaddata, cEffectHeatCapac[x] density];


along with the definition of assembler, quads, quaddata, etc., one can simply use

damping = RegionReactionMatrix[Ω, cEffectHeatCapac[x] density]


in the Do-loop.

• Looks great! But nothing works. Time is spent less than a second, but the output is nonsense. Mar 23, 2019 at 15:39
• Apparently, I forgot to DeployBoundaryConditions once before the Do-loop. Thanks for pointing that out! Mar 23, 2019 at 16:06
• Thanks, now it works. We can still speed up if we use option CompilationTarget -> "WVM" instead of CompilationTarget -> "C" (on my computer). How to get the dependence of temperature on time? Mar 23, 2019 at 16:49
• @Alex "We can still speed up if we use option CompilationTarget -> "WVM" instead of CompilationTarget -> "C" (on my computer)." That's odd and I don't know what to say about it. "How to get the dependence of temperature on time?" One has just to write x into successive rows of a matrix. OP did not do that in their post, so I presumed they were not interested in it. Or do you mean something else? Mar 23, 2019 at 17:10
• @OleksiiSemenov You're welcome. And yes, it is too bad that this method requires knowledge on how the actual matrices are assembled. I am pretty sure that the FEM developers would integrate such tools into NDSolve if they were given enough time and ressources. Alas, nonlinear PDEs in general can be very complicated and it isvery hard to write a robust user interface on the same level of abstraction as NDSolve. Mar 25, 2019 at 12:09

I managed to reduce the time by 2.5 times + I added the ability to display the temperature depending on the time. I used Do[] instead of For[] and Interpolation[] instead of Module[]. We can still speed up the code.

Needs["NDSolveFEM"];
Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
Lr = 2*10^-2;(*dimension of computational domain in r-direction*)Lz =
10^-2;(*dimension of computational domain in z-direction*)mesh =
ToElementMesh[FullRegion, {{0, Lr}, {0, Lz}},
MaxCellMeasure -> {"Length" -> Lr/50}, "MeshOrder" -> 1]
mesh["Wireframe"]
lambda = 22;(*heat conductivity*)density = 7200;(*density*)Cs = \
700;(*specific heat capacity of solid*)Cl = 780;(*specific heat \
capacity of liquid*)LatHeat =
272*10^3;(*latent heat of fusion*)Tliq = 1812;(*melting \
temperature*)MeltRange = 100;(*melting range*)To = 300;(*initial \
temperature*)SPow = 1000;(*source power*)R =
Lr/4;(*radius of heat source spot*)a = Log/R^2;
qo = (SPow*a)/Pi;
q[r_] := qo*Exp[-r^2*a];(*heat flux distribution*)tau =
10^-3;(*time step size*)ProcDur = 0.2;(*process duration*)
Heviside[x_, delta_] :=
Module[{res},
res = Piecewise[{{0,
Abs[x] < -delta}, {0.5*(1 + x/delta + 1/Pi*Sin[(Pi*x)/delta]),
Abs[x] <= delta}, {1, x > delta}}];
res]
HevisideDeriv[x_, delta_] :=
Module[{res},
res = Piecewise[{{0,
Abs[x] > delta}, {1/(2*delta)*(1 + Cos[(Pi*x)/delta]),
Abs[x] <= delta}}];
res]
EffectHeatCapac[tempr_] :=
Module[{phase}, phase = Heviside[tempr - Tliq, MeltRange/2];
Cs*(1 - phase) + Cl*phase +
LatHeat*HevisideDeriv[tempr - Tliq, 0.5*MeltRange]]
ehc = Interpolation[
Table[{x, EffectHeatCapac[x]}, {x, To - 100, 4000, 1}]];
ts = AbsoluteTime[];

NumTimeStep = Floor[ProcDur/tau];

vd = NDSolveVariableData[{"DependentVariables" -> {u},
"Space" -> {r, z}, "Time" -> t}];
sd = NDSolveSolutionData[{"Space",
"Time"} -> {ToNumericalRegion[mesh], 0.}];

DirichCond = DirichletCondition[u[t, r, z] == To, z == 0];
NeumCond = NeumannValue[q[r], z == Lz];
initBCs =
InitializeBoundaryConditions[vd, sd, {{DirichCond, NeumCond}}];
methodData = InitializePDEMethodData[vd, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];
xlast = Table[{To}, {methodData["DegreesOfFreedom"]}];
TemprField = ElementMeshInterpolation[{mesh}, xlast];
Do[(*(*Setting of PDE coefficients for linear \
problem*)pdeCoefficients=InitializePDECoefficients[vd,sd,\
"ConvectionCoefficients"\[Rule]{{{{-lambda/r,0}}}},\
"DiffusionCoefficients"\[Rule]{{-lambda*IdentityMatrix}},\
"DampingCoefficients"\[Rule]{{Cs*density}}];*)(*Setting of PDE \
coefficients for nonlinear problem*)
pdeCoefficients =
InitializePDECoefficients[vd, sd,
"ConvectionCoefficients" -> {{{{-(lambda/r), 0}}}},
"DiffusionCoefficients" -> {{-lambda*IdentityMatrix}},
"DampingCoefficients" -> {{ehc[TemprField[i - 1][r, z]]*density}}];
discretePDE = DiscretizePDE[pdeCoefficients, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];
A = damping/tau + stiffness;
x = LinearSolve[A, b,
Method -> {"Krylov", Method -> "BiCGSTAB",
"Preconditioner" -> "ILU0",
"StartingVector" -> Flatten[xlast, 1]}];
TemprField[i] = ElementMeshInterpolation[{mesh}, x];
xlast = x;, {i, 1, NumTimeStep}]
te = AbsoluteTime[];
te - ts


Here is the time for the old and new code 39.4973561 and 15.4960282 respectively (on my ASUS ZenBook).To further reduce the time, use the option MeshRefinementFunction:

f = Function[{vertices, area},
Block[{r, z}, {r, z} = Mean[vertices];
If[r^2 + (z - Lz)^2 <= (Lr/4)^2, area > (Lr/50)^2,
area > (Lr/
15)^2]]];
mesh =
ToElementMesh[FullRegion, {{0, Lr}, {0, Lz}}, "MeshOrder" -> 1,
MeshRefinementFunction -> f]
mesh["Wireframe"]


For this option time is 8.8878213

{ContourPlot[TemprField[NumTimeStep][r, z], {r, 0, Lr}, {z, 0, Lz},
PlotRange -> All, ColorFunction -> "TemperatureMap",
PlotLegends -> Automatic, FrameLabel -> Automatic],
ListPlot[Table[{tau*i, TemprField[i][.001, Lz]}, {i, 0,
NumTimeStep}], AxesLabel -> {"t", "T"}]} Thanks to Henrik Schumacher, we can still speed up the code. I slightly edited his code in case of using the "WVM" and to display the temperature field at each step.

Needs["NDSolveFEM"];
Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];

SetAttributes[AssemblyFunction, HoldAll];

Assembly::expected =
"Values list has 2 elements. Expected are 1 elements. Returning \
prototype.";

Assemble[pat_?MatrixQ, dims_, background_: 0.] :=
Module[{pa, c, ci, rp, pos},
pa = SparseArraySparseArraySort@SparseArray[pat -> _, dims];
rp = pa["RowPointers"];
ci = pa["ColumnIndices"];
c = Length[ci];
pos = cLookupAssemblyPositions[Range[c], rp, Flatten[ci], pat];
Module[{a},
a = <|"Dimensions" -> dims, "Positions" -> pos,
"RowPointers" -> rp, "ColumnIndices" -> ci,
"Background" -> background, "Length" -> c|>;
AssemblyFunction @@ {a}]];

AssemblyFunction /: a_AssemblyFunction[vals0_] :=
Module[{len, expected, dims, u, vals, dat}, dat = a[];
If[VectorQ[vals0], vals = vals0, vals = Flatten[vals0]];
len = Length[vals];
expected = Length[dat[["Positions"]]];
dims = dat[["Dimensions"]];
If[len === expected,
If[Length[dims] == 1, u = ConstantArray[0., dims[]];
u[[dat[["ColumnIndices"]]]] =
AssembleDenseVector[dat[["Positions"]], vals, {dat[["Length"]]}];
u, SparseArray @@ {Automatic, dims,
dat[["Background"]], {1, {dat[["RowPointers"]],
dat[["ColumnIndices"]]},
AssembleDenseVector[dat[["Positions"]],
vals, {dat[["Length"]]}]}}],
Message[Assembly::expected, expected, len];
Abort[]]];

cLookupAssemblyPositions =
Compile[{{vals, _Integer, 1}, {rp, _Integer, 1}, {ci, _Integer,
1}, {pat, _Integer, 1}},
Block[{k, c, i, j}, i = CompileGetElement[pat, 1];
j = CompileGetElement[pat, 2];
k = CompileGetElement[rp, i] + 1;
c = CompileGetElement[rp, i + 1];
While[k < c + 1 && CompileGetElement[ci, k] != j, ++k];
CompileGetElement[vals, k]], RuntimeAttributes -> {Listable},
Parallelization -> True, CompilationTarget -> "WVM",
RuntimeOptions -> "Speed"];

AssembleDenseVector =
Compile[{{ilist, _Integer, 1}, {values, _Real, 1}, {dims, _Integer,
1}}, Block[{A}, A = Table[0., {CompileGetElement[dims, 1]}];
Do[A[[CompileGetElement[ilist, i]]] +=
CompileGetElement[values, i], {i, 1, Length[values]}];
A], CompilationTarget -> "WVM", RuntimeOptions -> "Speed"];
Block[{q},
With[{code =
Flatten[Table[
Table[{CompileGetElement[q, i],
CompileGetElement[q, j]}, {i, 1, 4}], {j, 1, 4}], 1]},
Compile[{{q, _Integer, 1}}, code, CompilationTarget -> "WVM",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"]]];

Block[{dim, PP, UU, FF, p, u, f, integrant, x, \[Omega], localmatrix},
dim = 2;
PP = Table[CompileGetElement[P, i, j], {i, 1, 4}, {j, 1, dim}];
UU = Table[CompileGetElement[U, i], {i, 1, 4}];
FF = Table[CompileGetElement[F, i], {i, 1, 4}];
p = {s, t} \[Function] (PP[] (1 - s) + s PP[]) (1 - t) +
t (PP[] (1 - s) + s PP[]);
u = {s, t} \[Function] (UU[] (1 - s) + s UU[]) (1 - t) +
t (UU[] (1 - s) + s UU[]);
f = {s, t} \[Function] (FF[] (1 - s) + s FF[]) (1 - t) +
t (FF[] (1 - s) + s FF[]);
integrant = {s, t} \[Function]
Evaluate[f[s, t] u[s, t]^2 Abs[Det[D[p[s, t], {{s, t}, 1}]]]];
{x, \[Omega]} = Most[NIntegrateGaussRuleData[2, MachinePrecision]];
localmatrix =
1/2 D[Flatten[KroneckerProduct[\[Omega], \[Omega]]].integrant @@@
Tuples[x, 2], {UU, 2}];
getLocalDampingMatrices =
With[{code = localmatrix},
Compile[{{P, _Real, 2}, {F, _Real, 1}}, code,
CompilationTarget -> "WVM", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"]];];

Module[{fvals, fdata, localmatrices},
fvals = cEffectHeatCapac[Flatten[vals]]*density;
assembler[Flatten[localmatrices]]];
Lr = 2*10^-2;(*dimension of computational domain in r-direction*)Lz =
10^-2;(*dimension of computational domain in z-direction*)mesh =
ToElementMesh[FullRegion, {{0, Lr}, {0, Lz}},
MaxCellMeasure -> {"Length" -> Lr/50}, "MeshOrder" -> 1]
mesh["Wireframe"]

lambda = 22.;(*heat conductivity*)density = 7200.;(*density*)Cs = \
700.;(*specific heat capacity of solid*)Cl = 780.;(*specific heat \
capacity of liquid*)LatHeat =
272.*10^3;(*latent heat of fusion*)Tliq = 1812.;(*melting \
temperature*)MeltRange = 100.;(*melting range*)To = 300.;(*initial \
temperature*)SPow = 1000.;(*source power*)R =
Lr/4.;(*radius of heat source spot*)a = Log[100.]/R^2;
qo = (SPow*a)/Pi;
q[r_] := qo*Exp[-r^2*a];(*heat flux distribution*)tau =
10^-3;(*time step size*)ProcDur = 0.2;(*process duration*)
Heviside[x_, delta_] :=
Piecewise[{{0,
Abs[x] < -delta}, {0.5*(1 + x/delta + 1/Pi*Sin[(Pi*x)/delta]),
Abs[x] <= delta}, {1, x > delta}}];

HevisideDeriv[x_, delta_] :=
Piecewise[{{0,
Abs[x] > delta}, {1/(2*delta)*(1 + Cos[(Pi*x)/delta]),
Abs[x] <= delta}}];

EffectHeatCapac[tempr_] :=
Module[{phase}, phase = Heviside[tempr - Tliq, MeltRange/2];
Cs*(1 - phase) + Cl*phase +
LatHeat*HevisideDeriv[tempr - Tliq, 0.5*MeltRange]];
cHeaviside =
Compile[{{x, _Real}, {delta, _Real}},
Piecewise[{{0.,
Abs[x] < -delta}, {0.5*(1 + x/delta + 1./Pi*Sin[(Pi*x)/delta]),
Abs[x] <= delta}, {1., x > delta}}], CompilationTarget -> "WVM",
RuntimeAttributes -> {Listable}, Parallelization -> True];
cHeavisideDeriv =
Compile[{{x, _Real}, {delta, _Real}},
Piecewise[{{0.,
Abs[x] > delta}, {1./(2*delta)*(1. + Cos[(Pi*x)/delta]),
Abs[x] <= delta}}], CompilationTarget -> "WVM",
RuntimeAttributes -> {Listable}, Parallelization -> True];
cEffectHeatCapac[tempr_] :=
With[{phase = cHeaviside[tempr - Tliq, MeltRange/2]},
Cs*(1 - phase) + Cl*phase +
LatHeat*cHeavisideDeriv[tempr - Tliq, 0.5*MeltRange]];
u =.
vd = NDSolveVariableData[{"DependentVariables" -> {u},
"Space" -> {r, z}, "Time" -> t}];
sd = NDSolveSolutionData[{"Space",
"Time"} -> {ToNumericalRegion[mesh], 0.}];

DirichCond = DirichletCondition[u[t, r, z] == To, z == 0];
NeumCond = NeumannValue[q[r], z == Lz];
initBCs =
InitializeBoundaryConditions[vd, sd, {{DirichCond, NeumCond}}];
methodData = InitializePDEMethodData[vd, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

x0 = ConstantArray[To, {methodData["DegreesOfFreedom"]}];
TemprField = ElementMeshInterpolation[{mesh}, x0];
NumTimeStep = Floor[ProcDur/tau];

pdeCoefficients =
InitializePDECoefficients[vd, sd,
"ConvectionCoefficients" -> {{{{-(lambda/r), 0}}}},
"DiffusionCoefficients" -> {{-lambda*IdentityMatrix}},
"DampingCoefficients" -> {{EffectHeatCapac[
TemprField[r, z]] density}}];
discretePDE = DiscretizePDE[pdeCoefficients, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];
x = x0;
X = x;
taustiffness = tau stiffness;

assembler =
1], {1, 1} Length[mesh["Coordinates"]]];

A = damping + taustiffness;
x = LinearSolve[A, b,
Method -> {"Krylov", Method -> "BiCGSTAB",
"Preconditioner" -> "ILU0", "StartingVector" -> x,
"Tolerance" -> 0.00001}]; X[i] = x;, {i, 1, NumTimeStep}]; //
AbsoluteTiming // First


Here we have time 0.723424 and the temperature at each step

T[i_] := ElementMeshInterpolation[{mesh}, X[i]]

ContourPlot[T[NumTimeStep][r, z], {r, z} \[Element] mesh,
AspectRatio -> Lz/Lr, ColorFunction -> "TemperatureMap",
PlotLegends -> Automatic, PlotRange -> All, Contours -> 20]

ListPlot[Table[{i*tau, T[i][.001, Lz]}, {i, 0, NumTimeStep}]] • The section Transient PDE with Stationary Coefficients and Stationary Boundary Conditions shows how to set up a time dependent ElementMeshInterpolation Mar 25, 2019 at 7:06
• @user21 Here an explicit Euler is implemented for a non-linear equation. Standard code is implemented for a linear equation. Mar 25, 2019 at 11:25
• T[i_] := ElementMeshInterpolation[{mesh}, X[i]] is not the best way to do it. A better way is in the above mentioned tutorial. Mar 25, 2019 at 11:32
• @user21 Can you show with this example that standard code is faster? Mar 25, 2019 at 11:37
• It's not about faster, it's about better code. I am just trying to help you find better ways to code. Mar 25, 2019 at 11:43