Finite Element Method and ElementMeshInterpolation issue in v10.3

Question

Bug introduced in 10.3

---

I recently installed Mathematica v10.3 and checked a project I was working on few months ago under Mathematica v10.1.

In this project I solve a PDE with NDSolve and Finite Element Method over a 2D anulus region and I and ContourPlot that solution.

Under v10.1 (and v10.2) the solution appear like this:

Under v10.3 the solution appear like this:

The source code is exactly the same.

At the beginning I was convinced the problem was with NDSolve (and this question reflected that, I later edited). After more investigation I discovered the solution in the two releases is nearly the same.

You can reproduce the problem running:

Import["http://bit.ly/cp103issue"]

I searched on the documentation, news section, but I didn't find anything that could explain this or restore the previous behavior.

Any chance to understand the reason of the changed behavior and eventually restore the old?

UPDATE

The problem is not related to ContourPlot but to the InterpolatingFunction returned by NDSolve (thanks to JasonB). Probably the InterpolatingFunction is built under the cover by a call to ElementMeshInterpolation. Indeed the following code in Mathematica 10.3 returns result completely different from the results in Mathematica 10.2.

<<NDSolve`FEM`
{Ri,Ro}={3000000.,6371000.};
rgn=RegionDifference[Disk[{0.,0.},Ro],Disk[{0.,0.},Ri]];
mesh=ToElementMesh[rgn];
grid=mesh["Coordinates"];
vog=(Exp[-Norm[#]^2/Ri^2])&/@grid;
ListContourPlot[MapThread[Append,{grid,vog}],Epilog->{Point[grid]},PlotLegends->Automatic]
if=ElementMeshInterpolation[{mesh},vog];
ContourPlot[if[x,y],{x,y}\[Element]mesh,PlotLegends->Automatic,Epilog->{Point[grid]}]
TakeLargest[(if@@@grid-vog),20]

In 10.3:

In 10.2 (after som InterpolatingFunction::dmval message):

As you can see, in 10.2, the largest differences in "interpolation errors" are negligible, while in 10.3 they are even greather than the actual expected value.

The problem appear more specifically related with the "scale" of the mesh. Indeed if you change the sizes:

{Ri,Ro}={3000000.,6371000.}/10^6;

all apear fine:

Answer 1

Too long to make as a comment, so I'll post this non-answer as an answer.

I just want to point out that this definitely is not a ContourPlot problem, and it is absolutely a problem with the interpolating function.

So let's evaluate the following code in versions 10.2 and 10.31

Import["http://bit.ly/cp103issue"];
{sol1, sol2} = testSolution;
grid = testSolution[[1]]["Grid"];
vals1 = testSolution[[1]]["ValuesOnGrid"];
vals2 = testSolution[[2]]["ValuesOnGrid"];
norm = Norm /@ Transpose[{vals1, vals2}];
normgrid = Transpose[Append[Transpose[grid], norm]];

So I have extracted the grid points from the interpolating functions, as well as the values on those grid points. I have checked, and they are the same across both versions to within an accuracy of 10^-17. So we should be dealing with the same interpolating functions. But we aren't.

Here is the nail in the coffin,

ListPlot[sol1 @@@ grid - vals1]
ListPlot[sol2 @@@ grid - vals2]

Evaluating this code in version 10.2, we get

which looks great. In 10.3.1 we get

which is just plain horrible. So what we have here is some backsliding in the interpolating function resulting from an NDSolve call using a nonrectangular region.

Now here is something interesting, if you take those grid points and the values on the grid, and you make an interpolation function using Interpolation, then you get basically the same plot in both versions.

ContourPlot[
 Interpolation[Transpose[{grid, norm}], InterpolationOrder -> 1][x, 
  y], {x, y} ∈ testMesh]

Answer 2

Update: This is fixed in Version 11.0:

<< NDSolve`FEM`
{Ri, Ro} = {3000000., 6371000.};
rgn = RegionDifference[Disk[{0., 0.}, Ro], Disk[{0., 0.}, Ri]];
mesh = ToElementMesh[rgn];
grid = mesh["Coordinates"];
vog = (Exp[-Norm[#]^2/Ri^2]) & /@ grid;
(*ListContourPlot[MapThread[Append,{grid,vog}],Epilog\[Rule]{Point[\
grid]},PlotLegends\[Rule]Automatic]*)

if = ElementMeshInterpolation[{mesh}, vog];
(*ContourPlot[if[x,y],{x,y}\[Element]mesh,PlotLegends\[Rule]Automatic,\
Epilog\[Rule]{Point[grid]}]*)
TakeLargest[(if @@@ grid - vog), 20]

{4.996*10^-16, 4.85723*10^-16, 4.44089*10^-16, 3.88578*10^-16, 
 3.33067*10^-16, 3.33067*10^-16, 2.498*10^-16, 2.498*10^-16, 
 2.22045*10^-16, 2.22045*10^-16, 2.22045*10^-16, 2.22045*10^-16, 
 2.22045*10^-16, 1.94289*10^-16, 1.66533*10^-16, 1.66533*10^-16, 
 1.66533*10^-16, 1.66533*10^-16, 1.66533*10^-16, 1.66533*10^-16}

Old message:

This is a bug; you can workaround that by setting

SetSystemOptions[ "FiniteElementOptions" -> {"CacheInterpolationElements" -> False}]

Update:

It seems the default tolerance is not high enough for mesh over a large space compared to the default bounding box:

SetSystemOptions[
 "FiniteElementOptions" -> {"InterpolationToleranceFactor" -> 10^-8}]

This would allow you to still make use of cached elements which is faster for structured queries (which *Plot is not)

Sorry about that.