WhenEvent for several variables

Question

I want to solve a system of nonlinear two differential equations (say, in $\theta(x,t)$ and $\phi(x,t)$ with NDSolve, but I want to stop the evaluation when one of the functions ($\phi(x,t)$) reaches zero, or, for instance, a certain value $\bar{\phi}$.

Because then I would like to study how, where and (of course) when this pinch-off happens.

What is the best way of doing this, in your opinion?

Putting a constraint before running NDSolve?

I tried with WhenEvent, adding a line WhenEvent[\[Phi][x, t] == 0, Plot[\[Phi][x, T], {x, 0, L}] in the code, but it gives me

NDSolve::litarg:
To avoid possible ambiguity, the arguments of the dependent variable in  WhenEvent
(...)should literally match the independent variables

Here is the link to the shortest code I can manage to give:

https://www.dropbox.com/s/4t5o53lppb2itk8/mma3.nb

Answer 1

What Silvia wrote about the method of lines is basically true. The problem is that NDSolve tries as much as possible to "hide" the specifics about the discretization procedure so you can't use something like u[t,x] in a WhenEvent for a PDE as x would mean "any x" and that would be inifintely many tests. But for exactly that use-case there is the possibility to get a list of values at the actual discretization points instead which makes the formulation of the desired WhenEvent almost trivial and quite fast. The trick is to use "DiscretizedMonitorVariables" -> True for the "MethodOfLines" method. u[t,x] will then in fact be a vector of as many values as the method of lines uses in the spatial direction (and I think it should handle all cases with nontrivial spatial discretizations as well):

NDSolveValue[{
  D[u[t, x], t] == D[u[t, x], x, x],
  u[0, x] == 0,
  u[t, 0] == Sin[t],
  u[t, 5] == 0,
  WhenEvent[Min[u[t, x]] < -10^-4,
   tEnd = t; "StopIntegration"
   ]
  }, u, {t, 0, 10}, {x, 0, 5},
 Method -> {"MethodOfLines", "DiscretizedMonitorVariables" -> True}
 ]

For more information, see tutorial/NDSolveWhenEvents#24434 in the documentation center (or search for "Avoiding Wraparound in PDEs" here) and of course tutorial/NDSolveOverview is always a good place to learn about the details of NDSolve...

Answer 2

I haven't seen your particular PDEs, but I don't think your can numerically evaluate the target functions in the NDSolve scope, as the numerical solutions (expressed as InterpolatingFunction objects) hasn't been constructed.

If it were me, I may consider doing the plotting afterward:

sol = NDSolve[{
                y''[t] == -1, y[0] == 1, y'[0] == 0,
                WhenEvent[
                           y[t] == 0,
                           tEnd = t; "StopIntegration"
                         ]
              }, y, {t, 0, 2}];
Plot[y[t] /. sol, {t, 0, tEnd}]

Update:

I think I overlooked the question, and the real difficulty here might be related to the numerical method of lines, which is the way how NDSolve deals with this kind of PDEs.

According to the documentation, with the NML, the directions of boundary-value problem are considered as space dimensions, and the direction of initial-value problem is considerd as time dimension. All space dimensions are discretized, but the numerical integrations are only taken along the time dimension. So it's many lines start from some points in the space space at time zero, and evolve over time, thus the internal integrations are all one-dimensional and the "StopIntegration" action can only be taken along time. This is the reason you got the error in your comment.

One possible (and only partial) solution is specifying all space variables in WhenEvent, leaving only the time variable. Here is a simple demonstration:

Module[{sol},
       sol = NDSolve[{
                      D[u[t, x], t] == D[u[t, x], x, x],
                      u[0, x] == 0,
                      u[t, 0] == Sin[t], u[t, 5] == 0,

                      (* the space variable x, i.e. #, is fixed value: *)
                      WhenEvent[u[t, #] - .1 == 0,
                                tEnd = t; "StopIntegration"
                               ]

                      },
                      u, {t, 0, 10}, {x, 0, 5}];
       Plot3D[
              Evaluate[-(#/10) + u[t, x] /. sol], {t, 0, tEnd}, {x, 0, 5},
              PlotRange -> All, 
              PlotStyle -> ColorData["Rainbow"][1 - Rescale[#, {0, 5}]], 
              Mesh -> False]
      ] & /@ Range[0, 5, .5] // 
 Show[#, AxesLabel -> (Style[#, 20, Italic, Bold] & /@ {t, x, u}), Lighting -> "Neutral"] &