# Numerical instabilities of a convection-(non-)diffusion equation when shrinking from a square to a triangular domain

I am trying to evaluate a parameter-dependent indefinite integral using a PDE-based scheme I described here, and I'm having some trouble when I try and cut down the domain from a square to a triangle.

To be a bit more specific, suppose I'm trying to solve the differential equation $$\frac{\partial F}{\partial t}(t,t')=f(t,t')=\cos(t)-\cos(t')$$ under $F(0,t)=0$ for the triangular domain $t,t'\geq 0$, $t+t'\leq 50$. If I try to solve the problem on a square,

sol = F /. First[NDSolve[{
D[F[t, tt], t] == Cos[t] - Cos[tt]
, F[0, tt] == 0
}, f
, {t, 0, 50}, {tt, 0, 50}
]]


it evaluates fine, and plotting the function on the triangle returns

Plot3D[
sol[t, tt]
, {t, 0, 50}, {tt, 0, 50}
, AxesLabel -> {"t", "tt", ""}
, ImageSize -> 600
, PlotPoints -> 50
, RegionFunction -> Function[{t, tt, f}, t + tt < 50]
]


which is what it needs to return.

However, the problem starts if I try to solve directly in the triangular region, as

sol = F /. First[NDSolve[{
D[F[t, tt], t] == Cos[t] - Cos[tt]
, F[0, tt] == 0
}, f
, Element[{t, tt}, Triangle[{{0, 0}, {50, 0}, {0, 50}}]]
]]


This evaluates, returns an ugly function that I can plot,

and it returns the error message.

NDSolve::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

As I wrote here, this seems to be an indication that Mathematica thinks that the equation is a convection-diffusion equation with no diffusion in it (which, as long as it evaluates correctly, is just fine by me), but for the triangular region it suddenly somehow decides that there's no longer enough diffusion going on for things to work.

In other uses of this method, I've found it helpful to follow the error message's advice and add some artificial diffusion, in the form of a small multiple of the Laplacian on the right-hand side of the differential equation:

sol = F /. First[NDSolve[{
D[F[t, tt], t] == Cos[t] - Cos[tt] + c Laplacian[F[t, tt], {t, tt}]
, F[0, tt] == 0
}, f
, Element[{t, tt}, Triangle[{{0, 0}, {50, 0}, {0, 50}}]]
]]


The problem here, though, is that it requires extremely large values of the multiplier c (on the order of unity), which does change the behaviour of the solution way past acceptable bounds. Making c small, on the other hand, causes the error message

NDSolve::femcsp: The computed Peclet number is 6.841504974954965 and is larger than the mesh order (2), and the result may not be stable. Adding artificial diffusion may help.

which again sort of hints at Mathematica thinking that there's not enough diffusion going on, with the added bonus that there is something like mesh problems going on, which makes the change of behaviour on a smaller domain make more sense.

So, the question is, is there any way out of this conundrum? Mathematica was perfectly fine with the equation on a square domain; is there some way to re-embrace it on the triangular subset?

In particular, is there some way to make NDSolve choose a method that will work better for this setting? The method of lines looks tailor-made for this problem (i.e. thinking of the equation as a over-simple transport equation rather than a truncated convection-diffusion equation) but I've been unable to make it solve the problem.

One particular thing to note is that in general I will have relatively limited knowledge of the right-hand side integrand $f$, and it will in practice not be available for parameters outside of the triangular domain, so I do care about imposing the domain condition rather than simply computing the indefinite integral for a larger domain than I really need.

I only saw this question now, as it did not have the finite element tag. I hope an answer is still of interest.

Here is what happens: Since NDSolve gets a mesh for the entire region it can not deal with this problem as a time dependent problem. Most cases when I see this message the first thing to try is to separate the time domain ({t, 0, tEnd}) and the spatial domain (Element[{x},mesh]) and possibly adding an option (Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement"}}}).

Now, in this case this does not seam possible since the spatial domain changes. One way to deal with this is to add artificial diffusion. In order to avoid a too large Peclet number the mesh needs to be very fine. If you get a message about a too large Peclet number you need to refine the mesh. This then works:

Needs["NDSolveFEM"]
c = 0.01;
mesh = ToElementMesh[Triangle[{{0, 0}, {50, 0}, {0, 50}}],
MaxCellMeasure -> 0.01];
sol = NDSolveValue[{D[F[t, tt], t] ==
Cos[t] - Cos[tt] + c Laplacian[F[t, tt], {t, tt}], F[0, tt] == 0},
F, Element[{t, tt}, mesh]];

Plot3D[sol[t, tt], {t, 0, 50}, {tt, 0, 50},
AxesLabel -> {"t", "tt", ""}, ImageSize -> 600]


You can find some more information about this in the Finite Element Method Usage tutorial in the section about 'Stabilization of Convection-Dominated Equations'.