Where is the numerical solving breaking down?

Question

I am working with a set of three coupled reaction-diffusion PDEs, and for some parameter values I am getting some not so great solutions. I have been searching documentation and tutorials, and I have been trying to implement various methods and change various settings, but I feel like I am just blindly trying things without fully understanding what is going on. The equations might look intimidating, but I have supplied my code as well to help make everything reproducible.

As a summary, I was asking about why a bad solution is obtained for a certain parameter set. As I went through testing my problem to make this post I came across a second issue where the solution behavior also depended on the end time of integration, which seems odd to me. Why am I running into these issues?

---

Equations

The set of equations I am working with is as follows: $$\frac{\partial x}{\partial t}=\left(1+\frac{\gamma_{pL}\cdot y_T}{(1+\gamma_P\cdot x)^2}\right)^{-1}\cdot\left[\gamma_{D_P}\cdot\nabla^2x+\gamma_a\cdot x_I+\gamma_c\cdot\frac{\gamma_{pL}}{\gamma_P}\cdot\left(\frac{\gamma_P\cdot x}{1+\gamma_P\cdot x}\right)^2\cdot y_T\right]$$

$$\frac{\partial y_T}{\partial t}=-\gamma_c\cdot\frac{\gamma_P\cdot x}{1+\gamma_P\cdot x}\cdot y_T$$

$$\frac{\partial z}{\partial t}=\left(1+\frac{\gamma_R}{(1+\gamma_L\cdot z)^2}\right)^{-1}\cdot\left[\gamma_{D_L}\cdot\nabla^2z+\gamma_c\cdot\frac{\gamma_P\cdot x}{1+\gamma_P\cdot x}\cdot y_T\right]$$

where $x$, $y_T$, and $z$ are functions of $r$ (polar coordinate) and $t$ (time), $x_I$ is a source term given by $$x_I=\frac{1}{2\cdot\gamma_{D_P}\cdot t+1}\cdot\exp\left(-\frac{r^2}{2\cdot(2\cdot\gamma_{D_P}\cdot t+1)}-\gamma_a\cdot t\right)$$ and all $\gamma_i$ terms are parameters.

Initial conditions and boundary conditions are as follows: $$x(r,0)=z(r,0)=0,\ y_T(r,0)=1$$ $$\left.\frac{\partial x}{\partial t}\right|_{r=0}=\left.\frac{\partial z}{\partial t}\right|_{r=0}=0$$ $$x(\infty,0)=z(\infty,0)=0$$

In words, we start off with no $x$ and $z$, and $y_T$ starts at a maximum value of $1$ across all space. The no-flux boundary condition at the origin is because the system has radial symmetry, and we want $x$ and $z$ to go to $0$ at $\infty$.

---

Code

To tweak the above equations to work nicely in Mathematica (I think?), we need to use a "large" value of $r$ in place of $\infty$ and we need to use a point close to $r=0$ rather than $r=0$ itself (due to the $1/r$ term in the Laplacian in polar coordinates).

dr = .001; (*"small" r since we cannot use the origin in polar coordinates*)
ρFar = 100.; (*r that is "infinity"*)
τmax = 6.38;(*time to solve out to*)

The equations above:

xI = 1/(2*γDP*t + 1)*Exp[-(r-dr)^2/(2*(2*γDP*t + 1)) - γa*t];

dx = D[x[r, t],t] == (1 + (γpL*yT[r, t])/(1 + γP*x[r,t])^2)^-1*(γDP*Laplacian[x[r, t], {r, θ}, "Polar"] + γa*xI + γc*γpL/γP*((γP*x[r, t])/(1 + γP*x[r, t]))^2*yT[r, t]);

dyT = D[yT[r, t], t] == (-γc*γP*x[r, t]*yT[r, t])/(1 + γP*x[r, t]);

dz = D[z[r, t],t] == (1 + γR/(1 + γL*z[r, t])^2)^-1*(γDL*Laplacian[z[r, t], {r, θ},"Polar"] + (γc*γP*x[r, t]*yT[r, t])/(1 + γP*x[r, t]));

Intial/Boundary Conditions:

initx = x[r, 0] == 0.;

bc1x = (D[x[r, t], r] == 0) /. r -> dr;

bc2x = x[ρFar, t] == 0.;

inityT = yT[r, 0] == 1.;

initz = z[r, 0] == 0;

bc1z = (D[z[r, t], r] == 0) /. r -> dr;

bc2z = z[ρFar, t] == 0.;

Putting it all together:

deqns = {dx, dyT, dz, initx, bc1x, bc2x, inityT, initz, bc1z, bc2z};

Problem

My problem arises for certain parameter values (not all). I have also found that using

MaxStepSize -> .1

for NDSolveValue seems to help with some instances of my problem, but not all. An example of a "good solution" can be obtained with the following parameters

paramsGood = {γDP -> 0.001222, γDL -> 122.2, γa -> 1., γc -> 100., γP -> 0.01, γpL -> 0.01, γL -> 100., γR -> 0.01};


zSoln = NDSolveValue[deqns /. paramsGood, z, {r, dr, ρFar}, {t,0., τmax}, MaxStepSize -> .1];

Plot[zSoln[r, 2.6], {r, dr, 15}, PlotRange -> {-.002, .002}]

However, the following bad parameter set

paramsBad = {γDP -> 0.1222, γDL -> 122.2, γa -> 100., γc -> 1., γP -> 0.01, γpL -> 100., γL -> 1000., γR -> 100.};

gives an issue in the solution for $z$ (I have picked the time $t=2.6$ to show the bad solution, but the solution doesn't seem good from the start).

zSoln = NDSolveValue[deqns /. paramsBad, z, {r, dr, ρFar}, {t,0., τmax}, MaxStepSize -> .1];
Plot[zSoln[r, 2.6], {r, dr, 15}, PlotRange -> {-.002, .002}]

None of the dependent variables should be negative.

I have tried modifying the MaxStepSize, StartingStepSize, and specifying a spatial discretization as suggested in the documentation:

Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 1000}}}

but I cannot seem to fix this issue. Sometimes these things fix the problem, but then for other parameter values the issues come right back (for example, I think making the MinPoints option above larger fixes for this parameter set, but not for all of the. Additionally the solution ends up taking a long time). Why is this happening for the numerical solver, and what is the most efficient way to determine a well-behaved solution Note that I need to be able to solve these equations fairly quickly.

Second Problem

If I just change the end time to $t=6.5$ the solution ends up behaving differently. There is a very small response in $z$, but it at least looks well-behaved?

zSoln = NDSolveValue[deqns /. paramsBad, z, {r, dr, ρFar}, {t,0., 6.5}, MaxStepSize -> .1];
Plot[zSoln[r, 2.6], {r, dr, 15}, PlotRange -> {-.002, .002}]

Also, this issue does not occur for the "good parameter set"; the solution does not show any notable difference in choice of the final integration time.

So this really makes me question how do I make sure I am getting a legitimate solution here? Which settings should I employ to be confident in my solutions for a wide range of parameters?

Answer 1

For this problem, you must specify a solution method. Since the system of equations is nonlinear and the equation for y does not contain derivatives with respect to r, we must choose "MethodOfLines". It solves all problems right away. The code contains Quiet, since function xI is used outside the machine number, for example, 1/2.71828^726.112is calculated.

dr = .001;(*"small" r since we cannot use the origin in polar \
coordinates*)ρFar = 100.;(*r that is "infinity"*)τmax = \
6.38;(*time to solve out to*)
xI = 1/(2*γDP*t + 1)*
   Exp[-(r - dr)^2/(2*(2*γDP*t + 1)) - γa*t];

dx = D[x[r, t], 
    t] == (1 + (γpL*
          yT[r, t])/(1 + γP*x[r, t])^2)^-1*(γDP*
       Laplacian[x[r, t], {r, θ}, "Polar"] + γa*
       xI + γc*γpL/γP*((γP*
            x[r, t])/(1 + γP*x[r, t]))^2*yT[r, t]);

dyT = D[yT[r, t], 
    t] == (-γc*γP*x[r, t]*
      yT[r, t])/(1 + γP*x[r, t]);

dz = D[z[r, t], 
    t] == (1 + γR/(1 + γL*z[r, t])^2)^-1*(γDL*
       Laplacian[z[r, t], {r, θ}, 
        "Polar"] + (γc*γP*x[r, t]*
         yT[r, t])/(1 + γP*x[r, t]));
initx = x[r, 0] == 0.;

bc1x = (D[x[r, t], r] == 0) /. r -> dr;

bc2x = x[ρFar, t] == 0.;

inityT = yT[r, 0] == 1.;

initz = z[r, 0] == 0;

bc1z = (D[z[r, t], r] == 0) /. r -> dr;

bc2z = z[ρFar, t] == 0.;
deqns = {dx, dyT, dz, initx, bc1x, bc2x, inityT, initz, bc1z, bc2z};
paramsGood = {γDP -> 0.001222, γDL -> 
    122.2, γa -> 1., γc -> 100., γP -> 
    0.01, γpL -> 0.01, γL -> 100., γR -> 0.01};


zSoln = NDSolveValue[deqns /. paramsGood, 
    z, {r, dr, ρFar}, {t, 0., τmax}, 
    Method -> {"IndexReduction" -> Automatic, 
      "EquationSimplification" -> "Residual", 
      "PDEDiscretization" -> {"MethodOfLines", 
        "SpatialDiscretization" -> {"TensorProductGrid", 
          "MinPoints" -> 137, "MaxPoints" -> 137, 
          "DifferenceOrder" -> 2}}}] // Quiet;



paramsBad = {γDP -> 0.1222, γDL -> 122.2, γa -> 
    100., γc -> 1., γP -> 0.01, γpL -> 
    100., γL -> 1000., γR -> 100.};
zSoln1 = NDSolveValue[deqns /. paramsBad, 
    z, {r, dr, ρFar}, {t, 0., τmax}, 
    Method -> {"IndexReduction" -> Automatic, 
      "EquationSimplification" -> "Residual", 
      "PDEDiscretization" -> {"MethodOfLines", 
        "SpatialDiscretization" -> {"TensorProductGrid", 
          "MinPoints" -> 137, "MaxPoints" -> 137, 
          "DifferenceOrder" -> 2}}}] // Quiet;


zSoln2 = NDSolveValue[deqns /. paramsBad, 
    z, {r, dr, ρFar}, {t, 0., 6.5}, 
    Method -> {"IndexReduction" -> Automatic, 
      "EquationSimplification" -> "Residual", 
      "PDEDiscretization" -> {"MethodOfLines", 
        "SpatialDiscretization" -> {"TensorProductGrid", 
          "MinPoints" -> 137, "MaxPoints" -> 137, 
          "DifferenceOrder" -> 2}}}] // Quiet;


{Plot3D[zSoln[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All,
   PlotLabel -> "paramsGood", AxesLabel -> Automatic, Mesh -> None, 
  ColorFunction -> "Rainbow"], 
 Plot3D[zSoln1[r, t], {r, dr, 15}, {t, 0, τmax}, 
  PlotRange -> All, PlotLabel -> "paramsBad", AxesLabel -> Automatic, 
  ColorFunction -> "Rainbow"], 
 Plot3D[zSoln2[r, t], {r, dr, 15}, {t, 0, 6.5}, PlotRange -> All, 
  PlotLabel -> "paramsBad", AxesLabel -> Automatic, 
  ColorFunction -> "Rainbow"]}

Consider the following example

paramsGood = {γDP -> 0.01222, γDL -> 122.2, γa ->
     0.01, γc -> 100., γP -> 1., γpL -> 
    1., γL -> 1000., γR -> 100.} 

In this case, for "MaxPoints" -> 137, the solution goes into region z<0 (Fig. 2, the upper row on the left). We increase the number of points to 1037 and to 2037 (the top row is the center and to the right). The solution converges. However, with an increase of "DifferenceOrder" ->4 , the solution looks different (bottom row).This is due to artificial viscosity, which decreases with increasing "DifferenceOrder". I can recommend a combination of options "MinPoints" -> 1037, "MaxPoints" -> 1037, "DifferenceOrder" -> 2 for such cases (with the number of points you can experiment).

Answer 2

Some experimentation shows that the NDSolve problem is associated exclusively with z with no feedback to x or yT, because dz and dyT are independent of z. It is, therefore, convenient first to compute x and yT, and only then to compute z. For paramsGood,

deqns = {dx, dyT, initx, bc1x, bc2x, inityT};
deqnz = {dz, initz, bc1z, bc2z};

{xSoln, yTSoln} = NDSolveValue[deqns /. paramsGood, {x, yT}, {r, dr, ρFar}, 
    {t, 0., τmax}, MaxStepSize -> .1];
zSoln = NDSolveValue[deqnz /. paramsGood /. {x[r, t] -> xSoln[r, t], yT[r, t] -> yTSoln[r, t]}, 
    z, {r, dr, ρFar}, {t, 0., τmax}, MaxStepSize -> .1];
Plot3D[xSoln[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All, AxesLabel -> {r, t, x}, 
    PlotPoints -> 50, LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
    ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
 Plot3D[yTSoln[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All, AxesLabel -> {r, t, y}, 
    PlotPoints -> 50, LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
    ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[zSoln[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All, AxesLabel -> {r, t, z}, 
    PlotPoints -> 50, LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
    ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]

Visibly, the three computed veriables are well behaved everywhere. A slice through z at t = 2.6 reproduces the first plot in the question, as expected. The same code applied to paramsBad (but with PlotPoints -> 250) does not do so well.

Note the spikes near the origin in the plot of z. (The z computation in this case also produces a few underflow warnings, which I believe are insignificant.) Several constant-r slices through the z plot for very small r yield.

Plot[Table[zSoln[0.01 n, t], {n, 1, 9, 2}], {t, 0, .04}, PlotRange -> All, AxesLabel -> {t, z}]

On the other hand, a plot corresponding to the second plot in the question looks fine, although it may not be.

Plot[{zSoln[r, 2.6]}, {r, dr, 15}, PlotRange -> All, AxesLabel -> {r, z}]

In any case, the solution for z obtained here for paramsBad is much better than that in the question, probably because separating the computations of z allows NDSolve to optimize its automatic time-step determination for z alone.

Finally, let us repeat the computation but with MaxStepSize -> .01. Plots for x and yT appear identical to those just given, with the z plot now seemingly without error.

But, it isn't.

Plot[Table[zSoln[0.01 n, t], {n, 1, 9, 2}], {t, 0, .04}, PlotRange -> All, AxesLabel -> {t, z}]

Better, but not perfect. Note that the slice for t = 2.6 changes a bit in this case.

To answer the specific questions posed by the OP, first the errors obtained for paramsBad are due to lack of resolution, largely corrected by reducing MaxStepSize to 0.01. Still smaller MaxStepSize undoubtedly would produce a still more accurate result, but the computation becomes quite slow for 0.001. I also experimented some with very small StartingStepSize, but obtained little improvement. Second, the apparent improvement from the second to the third plots in the question are due to different time-steps used by NDSolve in the two cases, but plotting z in 3D shows that significant errors persist in the computation. Therefore, little credence should be given to the third plot in the question.

Answer 3

Fully NDSolve-based Solution

Adjustion for spatial step size together with temporal step size helps. I've used parameters mentioned in the comment for testing:

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

paramsBad = {γDP -> 0.01222, γDL -> 122.2, γa -> 0.01, γc -> 
    100., γP -> 1., γpL -> 1., γL -> 1000., γR -> 100.};

difforder = 4; points = 200;
test = NDSolveValue[
   deqns /. paramsBad, {x, yT, z}, {r, dr, ρFar}, {t, 0, τmax}, 
   Method -> mol[points, difforder], MaxStepSize -> {Automatic, 0.01}];
Plot3D[test[[-1]][r, t], {r, eps, 15}, {t, 0, τmax}, PlotRange -> All]

Partly NDSolve-based Solution

Another possible solution is discretizing the system in $r$ direction all by ourselves, in this case NDSolve seems to do a better job on choosing temporal step size so we don't need to adjust MaxStepSize option. Also, the singularity at $r=0$ is removed i.e. we can use dr = 0 in this method.

Notice pdetoode is used for discretization in the following program.

Clear[dr]
ρFar = 100.;
τmax = 6.38;
xI = 1/(2 γDP t + 1) Exp[-(r - dr)^2/(2 (2 γDP t + 1)) - γa t];    
dx = D[x[r, t], t] == (1 + (γpL yT[r, t])/(1 + γP x[r, t])^2)^-1 (γDP 
       Laplacian[x[r, t], {r, θ}, "Polar"] + γa 
       xI + γc γpL/γP ((γP x[r, t])/(1 + γP x[r, t]))^2 yT[r, t]);    
dyT = D[yT[r, t], t] == (-γc γP x[r, t] yT[r, t])/(1 + γP x[r, t]);
dz = D[z[r, t], t] == (1 + γR/(1 + γL z[r, t])^2)^-1 (γDL 
       Laplacian[z[r, t], {r, θ}, "Polar"] + (γc γP x[r, t] yT[r, t])/(1 + γP x[r, t]));    
initx = x[r, 0] == 0.;
bc1x = (D[x[r, t], r] == 0) /. r -> dr;
bc2x = x[ρFar, t] == 0.;
inityT = yT[r, 0] == 1.;
initz = z[r, 0] == 0;
bc1z = (D[z[r, t], r] == 0) /. r -> dr;
bc2z = z[ρFar, t] == 0.;
removesingularity = eq \[Function] r # & /@ Expand[eq] // Expand;
{newdx, newdz} = removesingularity /@ {dx, dz};
domain = {dr, ρFar};
points = 100; difforder = 2; grid = Array[# &, points, domain];
ptoofunc = pdetoode[{x, yT, z}[r, t], t, grid, difforder];
del = #[[2 ;; -2]] &;
ode = {del /@ ptoofunc@{newdx, newdz}, ptoofunc@dyT};
odeic = ptoofunc@{initx, inityT, initz};
odebc = With[{sf = 1}, ptoofunc@diffbc[t, sf]@{bc1x, bc2x, bc1z, bc2z}];
var = Outer[#[#2] &, {x, yT, z}, grid];
paramsBad = {γDP -> 0.01222, γDL -> 122.2, γa -> 0.01, γc -> 
    100., γP -> 1., γpL -> 1., γL -> 1000., γR -> 100.};
eps = 0;
sollst = NDSolveValue[{ode, odeic, odebc} /. paramsBad /. dr -> eps, 
   var /. dr -> eps, {t, 0, τmax}];
{solx, soly, solz} = rebuild[#, grid /. dr -> eps, 2] & /@ sollst;
Plot3D[solz[r, t], {r, 0, 15}, {t, 0, τmax}, PlotRange -> All]

Remark

1. As one can see, only 100 grid points are used. Increasing points won't cause significant difference.

2. Though difforder = 2 is chosen, difforder = 4 can be used. You may need to add Method -> {"EquationSimplification" -> Solve} option to NDSolve though.

3. Using dr = 0.001 won't cause any significant difference.

Answer 4

I think a key issue is the inability to resolve your source term, $x_I$, relative to your domain. Three parameter sets have been given through the course of the discussion. If we assume that we will use the same spatial and temporal domain for each parameter set, then a Plot3D of the source term reveals the relative feature size of the pulse.

Plot3D[xI /. paramsGood, {t, 0, τmax}, {r, dr, ρFar}, 
 PlotRange -> All, PlotLabel -> "Source Term paramsGood", 
 AxesLabel -> Automatic, Mesh -> None, ColorFunction -> "DarkBands", 
 PlotPoints -> 50, MaxRecursion -> 4]
Plot3D[xI /. paramsBad, {t, 0, τmax}, {r, dr, ρFar}, 
 PlotRange -> All, PlotLabel -> "Source Term paramsBad", 
 AxesLabel -> Automatic, Mesh -> None, ColorFunction -> "DarkBands", 
 PlotPoints -> 50, MaxRecursion -> 4]
Plot3D[xI /. paramsBad2, {t, 0, τmax}, {r, dr, ρFar}, 
 PlotRange -> All, PlotLabel -> "Source Term paramsBad2", 
 AxesLabel -> Automatic, Mesh -> None, ColorFunction -> "DarkBands", 
 PlotPoints -> 50, MaxRecursion -> 4]

As you can see, the first bad parameter set operates on a tiny portion of the domain. When the ratio of largest relevant feature to the smallest relevant feature that you want to capture is large, discretization can be difficult. In this case, we can use a helper function to create a non-uniform mesh that is very fine near $r=0$ and expands 100,000 times larger at the end of the domain with 200 grid points.

meshGrowth[x0_, xf_, n_, ratio_] := Module[{k, fac, delta},
  k = Log[ratio]/(n - 1);
  fac = Exp[k];
  delta = (xf - x0)/Sum[fac^(i - 1), {i, 1, n - 1}];
  N[{x0}~Join~(x0 + 
      delta Rest@
        FoldList[(#1 + #2) &, 0, 
         PowerRange[fac^0, fac^(n - 3), fac]])~Join~{xf}]
  ]
grid = meshGrowth[0.001, 100, 200, 100000];
ListLogPlot@grid

We can now apply this non-uniform grid to NDSolve. We also note that the source pulse is very short lived for the first parameters set. We will set the StartingStepSize and MaxStepFraction to small values to resolve the time domain. If we apply this to @bbgodfrey's approach, we obtain the following:

paramsGood = {γDP -> 0.001222, γDL -> 
    122.2, γa -> 1., γc -> 100., γP -> 
    0.01, γpL -> 0.01, γL -> 100., γR -> 0.01};
paramsBad = {γDP -> 0.1222, γDL -> 122.2, γa -> 
    100., γc -> 1., γP -> 0.01, γpL -> 
    100., γL -> 1000., γR -> 100.};
paramsBad2 = {γDP -> 0.01222, γDL -> 122.2, γa ->
     0.01, γc -> 100, γP -> 1, γpL -> 
    1., γL -> 1000., γR -> 100.};
dr = .001;(*"small" r since we cannot use the origin in polar \
coordinates*)
ρFar = 100.;(*r that is "infinity"*)
τmax = 6.38;(*time to solve out to*)
opts = (Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "Coordinates" -> {grid}}});
{xSolnG, yTSolnG} = 
  NDSolveValue[
   deqns /. paramsGood, {x, yT}, {r, dr, ρFar}, {t, 
    0., τmax}, opts, StartingStepSize -> 0.0001, 
   MaxStepFraction -> 1/400];
zSolnG = NDSolveValue[
   deqnz /. paramsGood /. {x[r, t] -> xSolnG[r, t], 
     yT[r, t] -> yTSolnG[r, t]}, 
   z, {r, dr, ρFar}, {t, 0., τmax}, opts, 
   StartingStepSize -> 0.0001, MaxStepFraction -> 1/400];
Plot3D[xSolnG[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All,
  AxesLabel -> {r, t, x}, PlotPoints -> 50, 
 LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
 ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[yTSolnG[r, t], {r, dr, 15}, {t, 0, τmax}, 
 PlotRange -> All, AxesLabel -> {r, t, y}, PlotPoints -> 50, 
 LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
 ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[zSolnG[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All,
  AxesLabel -> {r, t, z}, PlotPoints -> 50, 
 LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
 ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]

{xSolnB, yTSolnB} = 
  NDSolveValue[
   deqns /. paramsBad, {x, yT}, {r, dr, ρFar}, {t, 
    0., τmax}, opts, StartingStepSize -> 0.0001, 
   MaxStepFraction -> 1/400];
zSolnB = NDSolveValue[
   deqnz /. paramsBad /. {x[r, t] -> xSolnB[r, t], 
     yT[r, t] -> yTSolnB[r, t]}, 
   z, {r, dr, ρFar}, {t, 0., τmax}, opts, 
   StartingStepSize -> 0.0001, MaxStepFraction -> 1/400];
Plot3D[xSolnB[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All,
  AxesLabel -> {r, t, x}, PlotPoints -> 50, 
 LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
 ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[yTSolnB[r, t], {r, dr, 15}, {t, 0, τmax}, 
 PlotRange -> All, AxesLabel -> {r, t, y}, PlotPoints -> 50, 
 LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
 ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[zSolnB[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All,
  AxesLabel -> {r, t, z}, PlotPoints -> 50, 
 LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
 ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]

{xSolnB2, yTSolnB2} = 
  NDSolveValue[
   deqns /. paramsBad2, {x, yT}, {r, dr, ρFar}, {t, 
    0., τmax}, opts, StartingStepSize -> 0.0001, 
   MaxStepFraction -> 1/400];
zSolnB2 = 
  NDSolveValue[
   deqnz /. paramsBad2 /. {x[r, t] -> xSolnB2[r, t], 
     yT[r, t] -> yTSolnB2[r, t]}, 
   z, {r, dr, ρFar}, {t, 0., τmax}, opts, 
   StartingStepSize -> 0.0001, MaxStepFraction -> 1/400];
Plot3D[xSolnB2[r, t], {r, dr, 15}, {t, 0, τmax}, 
 PlotRange -> All, AxesLabel -> {r, t, x}, PlotPoints -> 50, 
 LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
 ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[yTSolnB2[r, t], {r, dr, 15}, {t, 0, τmax}, 
 PlotRange -> All, AxesLabel -> {r, t, y}, PlotPoints -> 50, 
 LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
 ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[zSolnB2[r, t], {r, dr, 15}, {t, 0, τmax}, 
 PlotRange -> All, AxesLabel -> {r, t, z}, PlotPoints -> 50, 
 LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
 ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]

Plot[Table[zSolnG[0.01 n, t], {n, 1, 9, 2}], {t, 0, .04}, 
 PlotRange -> All, AxesLabel -> {t, z}, PlotLabel -> "Good"]
Plot[Table[zSolnB[0.01 n, t], {n, 1, 9, 2}], {t, 0, .04}, 
 PlotRange -> All, AxesLabel -> {t, z}, PlotLabel -> "Bad"]
Plot[Table[zSolnB2[0.01 n, t], {n, 1, 9, 2}], {t, 0, .04}, 
 PlotRange -> All, AxesLabel -> {t, z}, PlotLabel -> "Bad2"]

We can also apply this non-uniform mesh to you original set of coupled equations for the first bad (i.e., sharpest pulse) parameter set and it seems to still hold together.

deqns = {dx, dyT, dz, initx, bc1x, bc2x, inityT, initz, bc1z, bc2z};
zSolnorig = 
  NDSolveValue[deqns /. paramsBad, 
   z, {r, dr, ρFar}, {t, 0., τmax}, opts, 
   StartingStepSize -> 0.0001, MaxStepFraction -> 1/400];
Plot[zSolnorig[r, 2.6], {r, dr, 15}, PlotRange -> All]
Plot3D[zSolnorig[r, t], {r, dr, 15}, {t, 0, τmax}, 
 PlotRange -> All, AxesLabel -> {r, t, z}, PlotPoints -> 50, 
 LabelStyle -> {Bold, Black, 15}, ImageSize -> Large, 
 ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]

By resolving the time and spatial domains of the source pulse, we have substantially eliminated negative values.