# WhenEvent used in a PDE to output independent variable

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[] but then I'm not sure how to get the output from ss[]?

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

• Can't you solve p1[z, t] == 20 after NDSolve is done? Or is there some reason you want to detect it during integration? – Michael E2 Jun 21 '15 at 18:44
• @MichaelE2 Yeah I think that should work equally well for now. I want to eventually package in Module but I guess I could still do that by this approach too. Anyway, I was attempting to use Solve just now but it's not working out. Would you mind showing me the string of code to solve for p1[z,t]==20? – E3labs Jun 21 '15 at 19:15
• I added it to my answer. – Michael E2 Jun 21 '15 at 19:47

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.

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.
*)

• Thanks a lot-All you guys here are great! The WhenEvent approach answers my original question so I consider this answered. FYI though, your addendum isn't working out. Just running zt20 I'm getting Position specification Missing["NotFound"] in Extract[p1["Grid"], Missing["NotFound"]] is not applicable and the output is Extract[p1["Grid"], Missing["NotFound"]]. Only thing I see is when I copy it over the p in FirstPosition's p<20 is green and not blue as in the post. Perhaps I'm missing your point about the extra braces. – E3labs Jun 21 '15 at 20:28
• @e3labs Oops, copied the wrong code. See the update. I usually prefer solving for the functions {p1, p2} instead of the value expressions {p1[z, t], p2[z, t]}, as it makes certain operations easier, such as computing residuals from pde and in this case, extracting properties of the InterpolatingFunction. (The blue in the post is due to Stackexchange; I have green too, which indicates that p is a local variable used to name the pattern.) – Michael E2 Jun 21 '15 at 20:57
• Thank you so much for the help Michael – E3labs Jun 21 '15 at 21:55