WhenEvent used in a PDE to output independent variable

Question

I want to use WhenEvent with a PDE after it has reached steady state. I'm posting a system of 2 equations (my real system has 6) and I'm solving the advection-diffusion-reaction 2nd order PDE.

First I set up the PDE:

pde = {
D[p1[z, t], t] == Dif D[p1[z, t], z, z] - v D[p1[z, t], z] - kp1 p1[z, t], 
D[p2[z, t], t] == Dif D[p2[z, t], z, z] - v D[p2[z, t], z] - kp2 p2[z, t] + kp1 p1[z, t]};

ic = {p1[z, 0] == p1in, p2[z, 0] == p2in};

bcs = {
v p1in == (v p1[z, t] - Dif D[p1[z, t], z] /. z -> 0), (D[p1[z, t],z] /. z -> l) == 0, 
v p2in == (v p2[z, t] - Dif D[p2[z, t], z] /. z -> 0), (D[p2[z, t],z] /. z -> l) == 0
};

params = {Dif -> 0.01, v -> 2, l -> 15, cin1 -> 10, p1in -> 30, p2in -> 60, kp1 -> 0.1, kp2 -> 0.05};

soli =NDSolve[
{pde, ic, bcs} /. params, {p1[z, t], p2[z, t]}, {z, 0, 15}, {t, 0, 100}, Method -> {"MethodOfLines", "SpatialDiscretization" ->{"TensorProductGrid", "MinPoints" -> 750}}];

I want to use WhenEvent to output the distance (i.e. z-variable) it takes for p1[z,t] (or p2[z,t]) to reach a chosen value. Something along the lines of:

WhenEvent[{p1[z,t]==20}, Print[distance along z]

My approach:

First I tried adding WhenEvent into soli like this:

soli =NDSolve[{pde, ic, bcs, WhenEvent[{p1[z, t] == 20}, {Print[z]}]} /. 
params, {p1[z, t], p2[z, t]}, {z, 0, 15}, {t, 0, 100}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 750}}];

That gave me all sorts of errors. My question is how do I get it to print out the distance along z that it take for p1[z,t] to decay to 20?

Attempt: Then I was wondering if I should discretize the solution so that I am working in 2D space instead of 3D space because it might be easier to extract the value? Something like:

ss = Table[Plot[p1[z, t] /. soli, {z, 0, 15}, PlotRange -> All], {t, 0, 10, .25}];

And then choose the last plot (i.e. so I'm pretty sure steady state is reached) with ss[[40]] but then I'm not sure how to get the output from ss[[40]]?

Part of me think I want to put it all into a module? Thanks in advance.

Answer 1

WhenEvent[cond, act] works in time, i.e., an event happens only when a time step causes the condition cond to change from False to True, save the special cases such as f == c described in the documentation. Those subtleties aside, the main thing to understand in using WhenEvent in the method of lines is what is substituted for the dependent variables such as p1[z, t]. If z is the spatial coordinate and t = t0 is the current value of t (the time front), then an interpolating function of z representing p1[z, t0] is substituted. You can deal with this as a function of z in any way you please.

Here is the OP's example modified appropriately:

soli = NDSolve[{pde, ic, bcs,
     WhenEvent[
      First@FindMinimum[{p1[z, t], 0 <= z <= 15}, {z, 10}] == 20, 
      Print[Join[
        FindRoot[(foo = p1[z, t]) == 20, {z, 0, 15}], {HoldPattern[t] -> t}]]]} /. params,
   {p1, p2}, {z, 0, 15}, {t, 0, 100}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 750}}];
(*  {z->11.8467, HoldPattern[t]->4.05457}  *)

Here is what is substituted for p1[z, t] in WhenEvent:

foo
(*  InterpolatingFunction[{{0., 15.}}, <>][z]  *)

Plot[foo, {z, 0, 15}]

Either of these may be used in place of First@FindMinimum[{p1[z, t], 0 <= z <= 15}, {z, 10}] == 20:

NMinValue[{p1[z, t], 0 <= z <= 15}, z] == 20  (* quite slow *)
Min[p1[z, t] /. z -> "ValuesOnGrid"] == 20          (* fast but less accurate *)

Both of these return {z->15., HoldPattern[t]->4.05457}, the same time value but different z value as above. Theoretically just taking the minimum of the interpolated values ("ValuesOnGrid") should not always give the same answer as actually minimizing the function, but here it agrees with NMinValue and not with FindMinimum. Global minimization (NMinValue, NMinimize) at each step by NDSolve can use up a lot of time.

Addendum

I suggested in a comment that determining when p1[z, t] == 20 might be done after the solution has been computed, when appropriate. It may be faster, too. Here is a way to proceed if soli is a solution returned by NDSolve as above (with the extra braces).

First find the point in the computed values where p1 first goes below 20. (Edit note: it is important that the expressions solved for in NDSolve are the functions {p1, p2}, not {p1[z, t], p2[z, t]} as in the OP's setup. The code below used the OP's first result throughout.)

zt20 = Extract[
  p1["Grid"] /. First[soli],
  FirstPosition[p1["ValuesOnGrid"] /. First[soli], p_ /; p < 20]
  ]
p1 @@ % /. soli
(*
  {8.02, 4.46224}
  {19.9818}
*)

Then we can use this as the starting point to find the least z where p1[z, t] == 20.

{min, cp} = FindMinimum[
  {z, 0 <= z <= 15 && 0 <= t <= 100 && p1[z, t] == 20} /. First[soli],
  Evaluate@Transpose[{{z, t}, zt20}]]
p1[z, t] /. cp /. First[soli]
(*
  {8.00134, {z -> 8.00134, t -> 4.436}}
  20.
*)