Hi everyone,

I am trying to solve a real-valued 2D advection-diffusion equation on the unit square. My Mathematica version is 12.1.1.0 on Windows 10. All seems to work well until the specified velocity coefficient functions get too complicated and reach some kind of internal threshold, in which case upon execution of NDSolveValue[] the kernel beeps and shuts down immediately. I would like to know if anyone can suggest a way forward to correcting this problem. Method details follow.

First a mesh is constructed:

W = ImplicitRegion[True,{{x,0,1},{y,0,1}}];
mesh = ToElementMesh[W, MaxCellMeasure->mcell];

The PDE operator is quite conventional

operator = Derivative[1,0,0][u[t,x,y]] - Dmol Derivative[0,2,0][u[t,x,y]] - Dmol Derivative[0,0,2][u[t,x,y]] + vx[t,x,y] Derivative[0,1,0][u[t,x,y]] + vy[t,x,y] Derivative[0,0,1][u[t,x,y]];

The boundary conditions are simple (natural Neumann boundaries are assumed elsewhere):

boundaryconditions = DirichletCondition[u[t,x,y]==0, x==1];

The initial condition is a simple Gaussian blob:

uInitial = Exp[-({x, y} - centroid).Inverse[s2].({x, y} - centroid)/2];

The problem parameters mcell, Dmol, vx, vy, centroid, s2 are specified as inputs. Once the various problem parameters are set, the solution is generated by the following statement:

uSolution = 
  NDSolveValue[{operator == 0, boundaryconditions, u[0, x, y] == uInitial}, 
   u, {t, 0, 1}, {x, y} \[Element] mesh, AccuracyGoal -> 10, InterpolationOrder -> All];

For simple cases with steady, uniform and homogeneous velocity coefficient functions {vx, vy} and low mesh Peclet numbers NDSolveValue returns accurate solutions quite readily for meshes with up to 100k triangles, which is all my RAM can handle. So far so good.

The problems start occurring when I specify more interesting velocity coefficient functions. The velocities of most interest to me are generated by high-resolution time-dependent interpolations over a 501x501 grid of points in the unit square. These data points are generated by an underlying FD scheme of dimension Nfd x Nfd. If Nfd is too large then, somehow, NDSolveValue[] simply crashes the kernel. The limit I have found is Nfd = 140, beyond which crashes occur no matter how coarse my FEM mesh is (even down to only 4k triangles). The "Why the Beep?" command reports simply that the kernel has encountered an exception.

I would like to have higher Nfd as this guarantees a physically accurate velocity field, but advecting the Gaussian blob in such a field just crashes immediately (before the first time step is made).

Questions:

1. I recognize that my preferred vx and vy are very complicated, time-dependent and spatially heterogeneous functions (although continuous and non-singular), but I don't understand how the dimension of the FD problem that specifies them is relevant to generating a kernel crash in a separate function call. There must be another problem, potentially memory related, or perhaps related to sampling of the vx,vy functions onto the FE mesh.

2. Is there any understandable tool that can help me debug a kernel crash?

Here are some example parameter values and a simply velocity specification that provides no problems to NDSolveValue[]:

s2 = IdentityMatrix[2]/625.0;
pos = {0.8,0.5};
Dmol = 0.0025;
mcell = 0.0004;
vx[t_,x_,y_] = -0.5;
vy[t_,x_,y_] = 0;

Cheers, Mike

PS: Further resources (notebook, results etc) are available if people are interested.