Aha, after playing with NDSolve for 10 years, I find the answer of my first question in this site!
Short Answer
I don't think my initial condition has that kind of fault: it's piecewise but smooth, right?
Yes, but prior error estimation of NDSolve doesn't only care about the smoothness of $T(0,x)$ (tes[0, x]) in this case. What's concerned is the smoothness of $\frac{\partial ^2T(0, x)}{\partial x^2} + e^{-1/T(0,x)}$ (D[tes[0, x], x, x] + Exp[-1/tes[0, x]]) i.e. right hand side (RHS) of the PDE at initial time.
Using higher order polynomial to build an initial condition (i.c.) with continuous second order derivative will resolve the problem, for example:
icexpr = Piecewise[{{2 - 1000 (0.1 - x)^3, 0 <= x <= 0.1},
{2 - 1000 (-0.9 + x)^3, 0.9 <= x <= 1}}, 2]
NDSolve[{D[tes[t, x], t] == D[tes[t, x], x, x] + Exp[-1/tes[t, x]],
tes[t, 0] == 1,
tes[t, 1] == 1, tes[0, x] == icexpr},
{tes[t, x]}, {t, 0, 100}, {x, 0, 1}]
(* Free of warning *)
Long Answer
First of all, let me untangle the sample a bit:
{xL = 0, xR = 1, dx = 0.1, T0 = 2, Tb = 1};
ic[n_][x_] :=
Piecewise[{{-dx^-n (dx - x)^n (T0 - Tb) + T0, xL <= x <= xL + dx},
{-dx^-n (x - (xR - dx))^n (T0 - Tb) + T0, xR - dx <= x <= xR}},
T0];
sample[n_, opt : OptionsPattern[]] :=
With[{T = T[t, x]},
NDSolve[{D[T, t] == D[T, x, x],
T == ic[n][x] /. t -> 0,
T == Tb /. {{x -> xL}, {x -> xR}}},
T, {t, 0, 100}, {x, xL, xR}, opt]]
sample[2]
NDSolve::mxsst
The $e^{-1/T}$ term (Exp[-1/tes[t, x]] in the question) has been removed because it turns out to be irrelevant to the problem.
To understand why the mxsst pops up, we need to know how the prior error estimation for PDE is done by NDSolve. The topic is discussed mostly in (but not limited to!) Spatial Error Estimates section of the obscure tutorial The Numerical Method of Lines. As a precondition, we also need to know how ODE(s) is pre-processed by NDSolve. I suggest reading this answer to understand the topic better………… OK, OK, I know it's non-trivial to read through all these, so let me try my best to retell the related parts in a easy-to-follow way.
When solving a time-dependent PDE (to be precise, PDE with independent variable that can be viewed as time), NDSolve will transform the PDE to a standard PDE, whose left hand side (LHS) is 1st order derivative respects to time. (Of course, for the specific problem discussed here, we don't need any further transform, because $\frac{\partial T}{\partial t}=\frac{\partial ^2T}{\partial x^2}$ is already in standard form. )
And the prior error estimation is built on the RHS of the obtained standard PDE.
The whole estimation process is a bit involved so I'd like not to talk too much about the details. We just need to remember: the RHS ($\frac{\partial ^2T}{\partial x^2}$ in our case) at initial time (t == 0 in our case) will be evaluated based on finite difference method (FDM), but the naive FDM won't be able to approximate a discontinuous function perfectly, because of the Gibbs phenomenon:
num = 101;
grid = If[$VersionNumber < 9, Table[x, {x, xL, xR, (xR - xL)/(num - 1)}],
Array[# &, num, {xL, xR}]];
rhsfunc[func_, grid_] := NDSolve`FiniteDifferenceDerivative[2, grid, func /@ grid]
rhs1 = rhsfunc[ic[2], grid];
Plot[ic[2]''[x], {x, xL, xR}, PlotRange -> All]~Show~
ListPlot[{grid, rhs1}\[Transpose], PlotStyle -> {Red}]
The error caused by Gibbs phenomenon contributes quite a bit to the error estimation of NDSolve. For num = 101, the estimated error is:
grid2 = If[$VersionNumber < 9, Table[x, {x, xL, xR, (xR - xL)/(2 num - 2)}],
Array[# &, 2 num - 1, {xL, xR}]];
rhs2 = rhsfunc[ic[2], grid2][[;; ;; 2]];
diff = Module[{lst = {rhs1, rhs2}},
lst[[All, ;; 3]] = 0.; lst[[All, -3 ;;]] = 0.;
Subtract @@ lst];
If[$VersionNumber < 11,
ListPlot[{grid, diff}\[Transpose] /. {x_, y_} /; Abs@y > 1 :>
Tooltip[{x, y}, TraditionalForm@"x" == N@x], PlotRange -> All] /.
Tooltip[{_, pts_}, label_] :> {pts, Text[label, pts[[1, 1]] + {0.07, 0.5}]},
ListPlot[{grid, diff}\[Transpose] /. {x_, y_} /; Abs@y > 1 :>
Callout[{x, y}, TraditionalForm@"x" == N@x], PlotRange -> All]]
p = If[$VersionNumber < 9, Infinity, 2];
NDSolve`ScaledVectorNorm[p, {0.00010264848819015053`, 0.00010264848819015053`},
NDSolve][diff, Internal`MaxAbs[rhs2, ic[2] /@ grid]]
(* In v9 or later: 3808.44 *)
(* In v8 or earlier: 27061.1 *)
Remark
"Are you sure you've reproduced the error estimation of NDSolve?"
Yes!:
norms = Trace[
sample[2, Method -> {"MethodOfLines",
"SpatialDiscretization" ->
{"TensorProductGrid",
"MinPoints" -> num}}],
HoldPattern@NDSolve`ScaledVectorNorm[__][__],
TraceInternal -> True] // Flatten // Union;
norms[[3]] // ReleaseHold
(* In v9 or later: 3808.44 *)
(* In v8 or earlier: 27061.1 *)
For any error estimation larger than 1, NDSolve will consider it as too large and increase number of grid aiming at a better estimation, but as we all know, the peak height of Gibbs phenomenon won't be relieved by denser grid, thus the mxsst warning finally comes up.
So, the problem won't show up when higher order polynomial is used to build the i.c. e.g. ic[3][x], because ic[3]''[x] is still continuous and is not influenced by Gibbs phenomenon:
Plot[ic[3]''[x], {x, xL, xR}, PlotRange -> All]~Show~
ListPlot[{grid, rhsfunc[ic[3], grid]}\[Transpose], PlotStyle -> {Red}]
sample[3]
(* Free of warning. *)
