Needs["NDSolve`FEM`"];
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];

Lr = 2*10^-2; (*dimension of computational domain in r-direction*)
Lz = 10^-2;   (*dimension of computational domain in z-direction*)
mesh = ToElementMesh[FullRegion[2], {{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_] := 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]
                                 ]

ts = AbsoluteTime[];
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[2]}},
"DampingCoefficients" -> {{Cs*density}}];
   *)

(*Setting of PDE coefficients for nonlinear problem*)

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

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

 A = damping/tau + stiffness;
 b = load + damping.xlast/tau;

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