Piecewise causes ndnum warning in NDSolve

Question

Bug introduced in 5.0, persisting through 13.2.

---

I try to solve the heat transfer equation with boundary condition that depends on time:

r0 = 0.75 10^-3;(*Beam spot size, m*)
ω = π ν;
ν = 1; (*pulse repetition rate, Hz*)

c = 1710;(*Heat capacity, W/(m·K)*)
ρ = 879; (*Density, kg/m^3*)
λ = 0.111;(*Heat conductivity, W/(m·K)*)
T0 = 300;(*initial temperature,K*)
T1 = 1000 - T0; (*Hot state temperature, K*)
Rm = 3 10^-3;(*Sample Radius, m*)
zm = 2 10^-3 ;(*Sample thickness, m*)   

 eq = D[T[R, z, t], t] == λ/(
    c ρ) (1/(R + 10^-20) D[R D[T[R, z, t], R], R] + 
      D[T[R, z, t], {z, 2}]);

init1 = T[R, z, 0] == T0;
bc1 = D[T[R, z, t], {R, 1}] == 0 /. R -> 0;
bc2 = D[T[R, z, t], {R, 1}] == 0 /. R -> Rm;
bc3 = T[R, z, t] == 
    Piecewise[{{T0 + T1 (1 - Abs@Sin[ω t])^500, 
       0 <= R <= r0}, {T0, True}}] /. z -> 0;
bc4 = T[R, z, t] == T0 /. z -> zm;

sol = NDSolveValue[{eq, init1, bc1, bc2, bc3, bc4}, 
  T[R, z, t], {R, 0, Rm}, {z, 0, zm}, {t, 0, 10},
  AccuracyGoal -> 30, MaxStepFraction -> 0.05]

The piecewise function defines the temperature at the left side of the z-interval as a short square pulses coming with frequency ω to central part of the disk.

However, solver returns the error-messages:

NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. NDSolveValue::ndnum: Encountered non-numerical value for a derivative at t == 0.

What's wrong? I've tried with FEM package but it produce even more error-messages :)

Answer 1

Several issues here.

1. As mentioned by Feyre, the i.c. and b.c. are inconsistent. I'll set "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100} in response. For more information check this post.

2. ndnum is caused by a bug of NDSolve (yeah it's confirmed, I reported it WRI before), similar problem has been observed in this and this post. This can be resolved by adding Simplify`PWToUnitStep@ before Piecewise.

3. The modeling for laser pulse is improper. Ideally DiracDelta is a possible choice, but currently NDSolve can't handle it properly, so we need to use an approximate one, for example:

   dirac[r_, a_] = Sqrt[a/Pi] Exp[-a r^2];
   

   Periodity can be achieved with Mod:

   Plot[dirac[Mod[t, 1, -1/2], 100], {t, 0, 10}, PlotRange -> All]
   

   

Here's the solution:

(* Previous code is not modified so I'd like to omit them in this post. *)
dirac[r_, a_] = Sqrt[a/Pi] Exp[-a r^2];

bc3 = T[R, z, t] == Simplify`PWToUnitStep@
         Piecewise[{{T0 + T1*dirac[Mod[t, Pi/ω, -2^(-1)], 100], 
        0 <= R <= r0}, {T0, True}}] /. z -> 0;
bc4 = T[R, z, t] == T0 /. z -> zm;

mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
        "MinPoints" -> n, "DifferenceOrder" -> o}}
mol[tf : False | True, sf_: Automatic] := {"MethodOfLines",
  "DifferentiateBoundaryConditions" -> {tf, "ScaleFactor" -> sf}}

sol = NDSolveValue[{eq, init1, bc1, bc2, bc3, bc4}, 
    T, {R, 0, Rm}, {z, 0, zm}, {t, 0, 10}, 
    Method -> Union[mol[60, 2], mol[True, 100]]]; // AbsoluteTiming

Animate[Plot3D[sol[r, z, t], {r, 0, Rm}, {z, 0, zm}, PlotRange -> {0, 4000}, 
  PerformanceGoal -> "Quality"], {t, 0, 10}]

Remark

1. NDSolve spits out eerr warning, but I think it's not a big problem since the error is small. Use more grid points will probably suppress the error, you can have a try if you have time.

2. I think "FiniteElement" will do a better job on solving this problem.

3. I took away AccuracyGoal because of the reason mentioned here.