Issues With WhenEvent in PDEs

Question

I am trying to solve the PDE below:

xInit = 1;
xMax = 10;
wid = 0.01;
T = 1000000;
Cc = 0.5

tMax = 1;
PDE = {D[c[x, t], t] == D[c[x, t], {x, 2}] + f[x, t], 
   D[f[x, t], t] == -f[x, t]/T, 
   f[x, 0] == (-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, 
   c[x, 0] == 0, (D[c[x, t], x] /. {x -> 0}) == 0, 
   WhenEvent[c[x, t] - Cc == 0, f[x, t] -> 1]};

When I do so with:

PDESol = NDSolve[PDE, {c[x, t], f[x, t]}, {t, 0, tMax}, {x, 0, xMax}, 
   MaxStepFraction -> 0.001];

I get the following error:

NDSolve::nbnum1: The function value -0.5+InterpolatingFunction[{{0.,10.}},{5,7,0,{1001},{6},0,0,0,0,Automatic,{},{},False},{{0.,0.01,0.02,<<45>>,0.48,0.49,<<951>>}},{Developer`PackedArrayForm,{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,<<952>>},{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,<<951>>}},{Automatic}][x]==0 is not True or False when the arguments are <<1>>. >>

I understand that this has something to do with the method of lines and/or how Mathematica is solving this PDE, but I don't understand what is going wrong and how to fix my specific case.

Thanks!

Answer 1

Too long for a comment, but also an answer to how to use WhenEvent in the "MethodOfLines". WhenEvent may not be the tool to use here.

WhenEvent is a time event. It occurs at a time $t=t_0$ for all $x$. If you want an event to happen depending on the spatial location $x$ and time $t$, such as implied by the condition c[x, t] - Cc == 0, then you cannot do this with WhenEvent.

Here's how one might use WhenEvent in the OP's PDE. I doubt it is the solution the OP wants, but the point of presenting it is to show how to use WhenEvent in a PDE.

As pointed out in the linked answer, the instances of c[x, t] in the WhenEvent code are replaced by an InterpolatingFunction with the following properties:

1. The "Grid" of the InterpolatingFunction is the spatial grid being used in NDSolve.
2. The "ValuesOnGrid" of the InterpolatingFunction are the values of c[x, t0] at t0, where t0 equals the time of the event, and for x ranging over the spatial grid.

In an action such as f[x, t] -> val, the expression val needs to evaluate to a value for f[x, t0] for x ranging over the spatial grid. Consequently in the case at hand, val needs to be vector of the same length as the spatial grid.

One further complication is that WhenEvent, while having the attribute HoldAll, evaluates some of its arguments in preparation for integrating the differential equation. So it seemed to me that I had to create functions that would not evaluate until the InterpolatingFunction for c[x, t] was substituted. One can see in the plots below that this works. Again, it might not be what the OP wants; however, it might show the OP whether or not WhenEvent will address her or his needs.

maxval[if : InterpolatingFunction[___][x_]] := 
  Max[if /. x -> "ValuesOnGrid"];  (* could be Min[] or some other function *)
fnew[if : InterpolatingFunction[___][x_]] := 
  ConstantArray[1, Length[if /. x -> "Grid"]];
PDE = {D[c[x, t], t] == D[c[x, t], {x, 2}] + f[x, t], 
   D[f[x, t], t] == -f[x, t]/T, 
   f[x, 0] == (-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, 
   c[x, 0] == 0, (D[c[x, t], x] /. {x -> 0}) == 0, 
   WhenEvent[maxval[c[x, t]] == Cc, f[x, t] -> fnew[c[x, t]]]};

sol = NDSolve[PDE, {c[x, t], f[x, t]}, {t, 0, tMax}, {x, 0, xMax}, 
  MaxStepFraction -> 0.001]

NDSolve::bcart: Warning: an insufficient number of boundary conditions...

Plot3D[c[x, t] /. sol, {x, 0, 10}, {t, 0, 1}, PlotRange -> All, 
 AxesLabel -> Automatic]

We can see the event happened around t = 0.61:

Plot3D[f[x, t] /. sol, {x, 0, 10}, {t, 0, 1}, PlotRange -> All, 
 AxesLabel -> Automatic]

IIRC, WhenEvent in PDEs is not extensively documented. The example in the docs for WhenEvent does not depend on the spatial variable(s). For this reason, I thought this explanation would helpful for the site. (Maybe it deserves its own question?)

Answer 2

The PDF can be simplified by integrating the equation for $f (x,t)$ independently of $c (x,t)$, and then completing the definition of $f (x,t)$ on the line $c(x,t)=Cc$. The result is shown below.

xInit = 1;
xMax = 10;
wid = 0.01;
T = 1000000;
Cc = 0.5;

tMax = 1;

PDE = {D[c[x, t], t] == 
    D[c[x, t], {x, 2}] + 
     If[c[x, t] == Cc, 1, 
      Exp[-t/T]*(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi], c[x, 0] == 0, 
   Derivative[1, 0][c][0, t] == 0, c[xMax, t] == 0};
sol = NDSolveValue[PDE, c, {t, 0, tMax}, {x, 0, xMax}, 
   MaxStepFraction -> 0.001];

f[x_, t_] := 
 If[sol[x, t] == Cc, 1, 
  Exp[-t/T]*(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi]
{ContourPlot[sol[x, t], {x, 0, 4}, {t, 0, 1}, Contours -> 20, 
  PlotRange -> All, ColorFunction -> Hue, PlotLegends -> Automatic, 
  PlotLabel -> "c(t,x)", FrameLabel -> {"x", "t"}], 
 ContourPlot[f[x, t], {x, 0, 2}, {t, 0, 1}, Contours -> 20, 
  PlotRange -> All, PlotLegends -> Automatic, 
  ColorFunction -> "Rainbow", FrameLabel -> {"x", "t"}, 
  PlotLabel -> "f(x,t)"], 
 Plot[f[x, 0], {x, 0, 2}, AxesLabel -> {"x", "f(x,0)"}]}

In the statement of the author can not even see the line c[x,t]-Cc==0. I added code to solve this problem.

xInit = 1;
xMax = 10;
wid = 0.01;
T = 1000000;
Cc = 0.5;

tMax = 1;
PDE = {D[c[x, t], t] == D[c[x, t], {x, 2}] + f[x, t], 
   D[f[x, t], t] == -f[x, t]/T, 
   f[x, 0] == (-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, c[x, 0] == 0, 
   Derivative[1, 0][c][0, t] == 0, c[xMax, t] == 0, 
   WhenEvent[c[x, t] - Cc == 0, f[x, t] -> 1]};
sol = NDSolve[PDE, {c[x, t], f[x, t]}, {t, 0, tMax}, {x, 0, xMax}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 80, "MaxPoints" -> 100, 
       "DifferenceOrder" -> "Pseudospectral"}}] // Quiet
{ContourPlot[c[x, t] /. sol, {x, 0, 4}, {t, 0, tMax}, Contours -> 20, 
  PlotRange -> All, ColorFunction -> Hue, PlotLegends -> Automatic, 
  PlotLabel -> "c(t,x)", FrameLabel -> {"x", "t"}], 
 ContourPlot[f[x, t] /. sol, {x, 0, 2}, {t, 0, tMax}, Contours -> 20, 
  PlotRange -> All, PlotLegends -> Automatic, 
  ColorFunction -> "Rainbow", FrameLabel -> {"x", "t"}, 
  PlotLabel -> "f(x,t)"], 
 Plot[{(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, f[x, t] /. sol /. t -> 0,
    f[x, t] /. sol /. t -> tMax}, {x, 0, 2}, AxesLabel -> {"x", "f"}]}

Generally speaking the optionWhenEvent[c[x, t] - Cc == 0, f[x, t] -> 1] does not give the desired effect, so I wrote another code in which on the line $c (x, t) = Cc$, the function is defined as $t = t0 (x)$.Then the solution is determined. So that we can clearly see the effect I put T = 10.

xInit = 1;
xMax = 10;
wid = 0.01;
T = 10;
Cc = 0.5;
tMax = 1;
PDE1 = {D[c[x, t], t] == 
    D[c[x, t], {x, 2}] + 
     If[t >= t0[x], Exp[-(t - t0[x])/T], 
      Exp[-t/T]*(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi], c[x, 0] == 0, 
   Derivative[1, 0][c][0, t] == 0, c[xMax, t] == 0};
PDE = {D[c[x, t], t] == D[c[x, t], {x, 2}] + f[x, t], 
   D[f[x, t], t] == -f[x, t]/T, 
   f[x, 0] == (-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, c[x, 0] == 0, 
   Derivative[1, 0][c][0, t] == 0, c[xMax, t] == 0};
sol1 = NDSolveValue[PDE, c, {t, 0, tMax}, {x, 0, xMax}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 800, "MaxPoints" -> 1000, 
       "DifferenceOrder" -> "Pseudospectral"}}];

t0 = Interpolation[
  Table[{x, t /. FindRoot[sol1[x, t] == Cc, {t, .5}]}, {x, 0, 2, .1}]]
Sol2 = NDSolveValue[PDE1, c, {t, 0, tMax}, {x, 0, xMax}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 800, "MaxPoints" -> 1000, 
       "DifferenceOrder" -> "Pseudospectral"}}];

f1[x_, t_] := 
 If[t >= t0[x], Exp[-(t - t0[x])/T], 
  Exp[-t/T]*(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi]
{ContourPlot[Sol2[x, t], {x, 0, 4}, {t, 0, tMax}, Contours -> 20, 
  PlotRange -> All, ColorFunction -> Hue, PlotLegends -> Automatic, 
  PlotLabel -> "c(t,x)", FrameLabel -> {"x", "t"}], 
 ContourPlot[f1[x, t], {x, 0, 2}, {t, 0, tMax}, Contours -> 20, 
  PlotRange -> All, PlotLegends -> Automatic, 
  ColorFunction -> "Rainbow", FrameLabel -> {"x", "t"}, 
  PlotLabel -> "f(x,t)"], 
 Plot[{(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, f1[x, t] /. t -> 0, 
   f1[x, t] /. t -> tMax}, {x, 0, 2}, AxesLabel -> {"x", "f"}]}