Message Boards Message Boards


NDSolveValue: crashes kernel for FEM model

Posted 8 months ago
2 Replies
0 Total Likes

Hi everyone,

I am trying to solve a real-valued 2D advection-diffusion equation on the unit square. My Mathematica version is 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).


  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.

2 Replies

Can you post a non-crashing and a crashing example? I can run it in a newer version of Mathematica on a Mac and tell you if it crashes. I might also be able to give you some insight. Alternatively, you can send the examples to tech support -- they are always interested in kernel crashes.



Posted 7 months ago

Hi Neil, many thanks for your reply and your interest.

I was responding to you by developing a minimal example notebook to show (i) correct behaviour and (ii) the kernel crash. In doing so I extracted all the velocity solution method sections (which were very complicated) from the notebook and just read in the pre-calculated velocity solution from a dump file. Interestingly, as soon as I did this step, the kernel crash effect disappeared!

When dumping the solution I was very uncritical. After generating the velocity field I simply did the following:


to save all entities in the Global context. Then in my blob Notebook just before the advection-diffusion code I read in the prior Global context using


Then I go ahead to set parameter values and execute the NDSolveValue[] function as before. No problems!

As you can imagine I am both relieved and disturbed that this seems to avoid the kernel crash. Relieved because I can go ahead and calculate solutions with this work-around, but disturbed because I am no closer to understanding what the issue was in the first place.

Over the years I have encountered several issues where a notebook might get invisibly corrupted somehow, and I have had to manually copy and paste cells into a new notebook to get a clean working copy. Maybe my travails here are simply a different manifestation of this weirdness?

Cheers, Mike

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract