# Getting rid of spikes in the PDE solution

Bug introduced in 10.0 and fixed in 10.3

Note: In 10.0, Rationalize[fd, 0] was needed or mesh generation would fail.

Preamble: I am solving a PDE in a domain representing a rectangle 10X10 with several circular holes, say, like this one (n is the number and r is the radius of the holes):

Needs["NDSolveFEM"];
r = 1;
rct = Rectangle[{0, 0}, {10, 10}];
rnd[n_Integer] :=
Transpose[{SeedRandom[RandomInteger[1000]];
RandomReal[{r, 10 - r}, n], SeedRandom[RandomInteger[2000]];
RandomReal[{r, 10 - r}, n]}];

n = 10;
fd = Fold[RegionDifference, rct, Disk[#, r] & /@ rnd[n]];
mesh = ToElementMesh[fd];
Show[mesh["Wireframe"], ImageSize -> 200]


You should see this meshed domain:

The code above generates randomly placed holes, so playing with it one gets a number of domains with the same number and sizes of the holes. This can be used below for experimenting.

The PDE represents a plane Laplace equation \text{$\Delta$u}=0 with the Dirichlet boundary conditions u(x,0)=0 and u(x,10)=1.

In addition the solution along with its 3D plot is wrapped below by Manipulate. The latter enables one to play with the parameters controlling the mesh quality and size: MaxCellMeasure, mcm, and MeshQualityGoal, mqg:

    Clear[mesh];
Manipulate[
mesh = ToElementMesh[fd, "MaxCellMeasure" -> mcm,
"MeshQualityGoal" -> mqg];
slv = NDSolveValue[{D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] == 0,
DirichletCondition[u[x, y] == 1, y == 10],
DirichletCondition[u[x, y] == 0, y == 0]},
u[x, y], {x, y} \[Element] mesh];
Plot3D[Evaluate[slv], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow", PlotRange -> {-0.5, 1.5},
ClippingStyle -> Automatic,
AxesLabel -> {Style["x", 18, Italic], Style["y", 18, Italic],
Style["\[CurlyPhi]", 18, Italic]}, ImageSize -> 400],
{{mcm, 0.1}, Range[0.005, 0.2, 0.005]}, {{mqg, 1}, Range[0.2, 1, 0.2]}
]


After evaluation one should see this:

Now, one can see spikes in the solution. Since it is the Laplace equation, the inhomogeneities, if any, should have the lateral size comparable to the size of a typical geometric feature that has generated them. By inspecting the image one finds that it is not the case. At least, often not the case. Indeed, playing with the parameters mcm and mqg one can obtain a different spikes distribution. Here is what I see at mcm=0.05 and mqg=1:

It is, therefore, an erroneous result, at least, in most cases.

The further observation is that playing with mcm and mqg one can influence the number of spikes and, eventually, get rid of them, but I could not find any general rule leading to decreasing of the number of spikes this way, let alone removing them. Say, increasing or decreasing mqg may in some cases lead to the decreasing of the spikes number, while in the others to their increasing. The same effect gives the variation of mcm.

By further playing with the code above I observe that at small n values (n=1 to 3) there were no, or only a few spikes, but they increase in number and amplitude with the increase of n.

My future aim is to evaluate an integral of the squared gradient of u over the whole mesh. Here all contributions of spikes are summed up, and give sometimes unrealistic values. The whole calculation, therefore, appears to be unreliable.

My question: do you see the way to get rid of spikes after the solution is obtained?

Later edit: To address the question of Michael E2: "You say that the integral of the gradient squared sometimes gives unrealistic values. Can you give an example for us to explore? "

The expectation is that the values of the integral should be below unity. It is not excluded, however, that some of them might be slightly above unity. As it is shown below there are weird values much greater than unity.

This:

    ProgressIndicator[Dynamic[q], {1, 10}]
lstFull = Table[{n, Table[Module[{r, fd, mesh, slv},
q = n;
r = 1;
fd = Fold[RegionDifference, rct, Disk[#, r] & /@ rnd[n]];
mesh =
ToElementMesh[fd, "MaxCellMeasure" -> r/10,
"MeshQualityGoal" -> Automatic];
slv = NDSolveValue[{D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] == 0,
DirichletCondition[u[x, y] == 1, y == 10],
DirichletCondition[u[x, y] == 0, y == 0]},
u[x, y], {x, y} \[Element] mesh];
(NIntegrate[
Evaluate[D[slv, x]^2 + D[slv, y]^2], {x, y} \[Element] mesh])
], {100}]}, {n, 1, 10}];


makes a list (entitled "lstFull") with the structure {n, sublist} where the sublist contains 100 values of the integral in question obtained with 100 different configurations of n holes. For shortness I only varied n from 1 to 10. For the same sake I used the mesh "MaxCellMeasure" -> r/10, while I also tried the finer mesh and also tried the mesh order set to 2 with the same effect. The whole list lstFull is much too long to give it here, I will only select the weird terms:

g[{x_, y_}] := {x, Select[y, # > 100 &]};
Map[g, lstFull]


returning

(*  {{1, {}}, {2, {}}, {3, {}}, {4, {2.17472*10^6, 131807., 131807.,
131807., 131807., 131807.,
131807.}}, {5, {4.1142*10^6}}, {6, {}}, {7, {8.61214*10^11}}, {8, \
{3281.12, 129.226, 138.23, 138.23, 138.23, 138.23, 138.23, 138.23,
138.23}}, {9, {204.573, 414.913, 344448., 204.573, 414.913,
344448., 204.573, 414.913, 344448., 204.573}}, {10, {2023.6,
100.345, 172., 2.18092*10^6, 266.877, 63376.7, 2034.07, 2034.07,
2034.07, 2034.07, 2034.07, 2034.07, 2034.07, 2034.07, 2034.07}}}   *)


So, one can find the values as great as ~10^6 and ~10^11. I would say, however, that even the values as high as ~100 are unrealistic.

• They are interpolation or extrapolation problems. Head[slv]["ValuesOnGrid"] // SystemMinMax returns {0., 1.}. Most of the errors are on the boundary, but some are not which confuses me. If they were all on the boundary, one guess would be that Plot3D does not respect the quadratic mesh, but linearizes it and evaluates slv at points outside the mesh. (Playing with your example just crashed my Mma front end, so I thought I'd summarize where I thought I had gotten to.) – Michael E2 Jul 14 '15 at 15:29
• You say that the integral of the gradient squared sometimes gives unrealistic values. Can you give an example for us to explore? – Michael E2 Jul 14 '15 at 16:36
• your code in 'slv=' is not typeset properly? – chris Jul 15 '15 at 8:12
• @chris I do not see, where. Would be great, if you point out where precisely. – Alexei Boulbitch Jul 15 '15 at 8:21
• @Michael E2 Yes, it is clear that it is related to interpolation. My problem is how to make such that these errors do not interfere with the results. Is there, for example, a post processing possibility to get rid of spikes? For your second question please see the later edit. – Alexei Boulbitch Jul 15 '15 at 8:27

This is a bug (up until version 10.2) in the interpolation code when there are concave curved elements for unstructured query points.

Workaround: To work around it one can use the either a "MeshOrder"->1 or set the "ImproveBoundaryPosition"->False both of these get rid of the spikes.

Do[
r = 1;
rct = Rectangle[{0, 0}, {10, 10}];
rnd[n_Integer] :=
Transpose[{RandomReal[{r, 10 - r}, n],
RandomReal[{r, 10 - r}, n]}];

n = 10;
fd = Fold[RegionDifference, rct, Disk[#, r] & /@ rnd[n]];

mesh = ToElementMesh[fd
(*,"MeshOrder"\[Rule]1*)
, "ImproveBoundaryPosition" -> False
];
slv = NDSolveValue[{D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] == 0,
DirichletCondition[u[x, y] == 1, y == 10],
DirichletCondition[u[x, y] == 0, y == 0]},
u, {x, y} \[Element] mesh(*,
"ExtrapolationHandler" -> {(Indeterminate&),
"WarningMessage"\[Rule]False}*)];
Print[Plot3D[slv[x, y], {x, y} \[Element] mesh, PlotRange -> All]]
, {10}]


Lowering the mesh order will get a less accurate result, which could be compensated by using more elements. Not improving the boundary position will avoid inserting curved elements at the region boundary. This will result in a less accurate boundary approximation which can be compensated by refining the boundayr with "MaxBoundaryCellMeasure". My recomentation would be to go with "ImproveBoundaryPosition"->False if you are hit by this.

How does this happen: To understand what is going on one needs to understand how the interpolation algorithm works in principal. First a nearest algorithm is used to find a close element. That element is tested if the query point is in that element. If that is the case all is good and we interpolate. If it's not then the query of wheather a point is in the element or not will also indicate a direction in which to searach for another element that may contain the query point. This is done until the element if found or it is figured out that the query point is outside the region. The element which contained the last query point is cached. Since the likehoot that the element if queried again at a slightly different place is high. The check if a query point is in a concave curved element is complicated and has a fault if the query point is far from the actual element, as can happen if the last cached element is far from the next query point. That's why this is not noticed in sturtured queries and only during unstructured querying like done during plotting. If someone has ideas for an efficient and robust test if a point is inside a curved triangle, let me know in the comments. I am sorry for trouble this causes.

• Thank you, I am just in the process of checking the approach with the "ImproveBoundaryPosition" -> False  for different values of the holes radii and it presently looks pretty good: no spikes. – Alexei Boulbitch Jul 17 '15 at 9:30

There is the as yet undocumented FEM function (see this answer, or read the user-interface with Definition after clearing the ReadProtected attribute)

ElementMeshPlot3D[valuesOnMeshCoordinates, elementMesh, options]


We can use this to plot the surface (note Head[slv] is just the pure InterpolatingFunction of the solution so we can apply to methods; see 1, 2):

Needs["NDSolveFEM"];
Axes -> True, ColorFunction -> (ColorData[{"Rainbow", {-0.5, 1.5}}]),
PlotRange -> {-0.5, 1.5}, Options[Plot3D, BoxRatios]]


Ironically it does not do meshing of the plot.

Update: Checking the interpolation

This was confusing. Originally I tested slv on a regular grid

pts = Tuples[Range[0., 10., 0.01], 2]       (* input points *)


and everything looked perfectly good. But then I tried random points and got some strange results. Update -- I forgot to include this: I added an "ExtrapolationHandler" option so that the interpolation should return Indeterminate when an input point lies outside the mesh.

slv = NDSolveValue[{D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] == 0,
DirichletCondition[u[x, y] == 1, y == 10],
DirichletCondition[u[x, y] == 0, y == 0]},
u[x, y], {x, y} ∈ mesh,
"ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}]

SeedRandom[0];
pts = RandomReal[{0, 10}, {100000, 2}];   (* input points *)
vals = Head[slv] @@@ pts;                 (* solution values on input points *)
badpos = Position[vals, Indeterminate];   (* Indeterminate => points not in domain *)
goodpts = Delete[pts, badpos];            (* Delete the points not in domain *)
goodvals = DeveloperToPackedArray@
Delete[vals, badpos];                   (* Delete the corresponding values *)
(* debug info *)
foopts =              (* Pick the points with solutions values outside 0, 1 *)
Pick[goodpts, UnitStep[goodvals] UnitStep[1 - goodvals], 0];
foovals =             (* Save the corresponding original values of the solution *)
Pick[goodvals, UnitStep[goodvals] UnitStep[1 - goodvals], 0];
foorevals =           (* Reevaluate solution on points with excessive values *)
(* end debug *)
Graphics3D[        (* plots points vs. solution value, provided there was a value *)
GraphicsComplex[
Transpose[Append[Transpose@goodpts, goodvals]],
{PointSize[Small],
Point[Range@Length@goodpts,
VertexColors -> ColorData[{"Rainbow", {-0.5, 1.5}}] /@ goodvals]
}], Axes -> True, BoxRatios -> {1, 1, 1}
]
]


The values of the solution at some of the random points are excessively large in magnitude and lead to points lying well away from the surface. In some cases, the points are not even in the domain mesh.

The misbehaving points seem randomly placed with respect to the mesh:

Show[
mesh["Wireframe"],
Graphics[{PointSize[Medium],
Point[foopts,
VertexColors ->   (* color the points with excessive values according to reevaluation*)
(foorevals /.
{Indeterminate -> Red,  (* originally inside, reevaluated outside domain *)
_?NumericQ -> Blue})]  (* inside domain on both evaluations *)
}]]


Weirdest is that if slv is re-evaluated on the bad points (foorevals = Head[slv] @@@ foopts), I get different values, unless someone can point out a mistake (see also the graphics above):

foovals     (* original evaluation of solution was outside interval {0, 1} *)
foorevals   (* re-evaluation of solution is reasonable *)
(*
foovals yields (all numeric outside {0, 1}, no Indeterminate)
{6.20318, -1.15456, -0.225189, 4.75979, -15.0127, -0.655401, \
-0.17161, -1.10463, 2.08943, -2.66222, 26.6109, 41.465, -9.12591, \
-2.2055, 1.81575, 3.94393, 2.01316, 10.4552, 22.4335, 12.0272, \
59.6386, 23.2775, 5.53432, 7.1062, -3.05542, -22.1865, -2.22477, \
-40.5778, 2.79623, -0.294799, 41.8734, 19.2247, -3.87972, 1.61702, \
4.37806, 12.1792, 1.95929, -30.1799, -4.32195, 12.9452, 30.5671, \
3.05843, 36.1302, 9.3391, 11.4406, -41.6207, 9.02034, -1.36818, \
-2.11299, -6.87225, -57.8951, -32.7291, 4.83439, 2.0149, 14.4227, \
-5.46179, -14.3458, -0.422857, 8.94631, 2.80505, 1.19526, -0.483842, \
-0.51311, 2.12031, 8.81922, 1.89209, -63.8548, -20.6707, 1.71642, \
-0.0584718, 13.4679, 2.89107}

foorevals yields (either inside {0, 1} or outside domain = Indeterminate)
{0.861886, 0.00673629, 0.0370507, 0.78554, 0.843989, Indeterminate, \
Indeterminate, 0.00226378, 0.817398, 0.648706, 0.779712, \
Indeterminate, 0.833617, 0.644034, 0.829276, 0.702965, 0.989782, \
0.987569, 0.943198, 0.799872, 0.945101, 0.951939, 0.961857, 0.214184, \
0.699705, Indeterminate, 0.866424, 0.987908, 0.681963, 0.00692621, \
0.637715, Indeterminate, 0.0111853, 0.858184, 0.972909, 0.858585, \
0.659178, 0.93633, 0.752761, 0.573878, 0.358457, 0.732126, 0.947696, \
0.12764, 0.686666, 0.997436, 0.656519, 0.718724, 0.717775, 0.859168, \
0.965433, 0.870962, 0.985584, Indeterminate, 0.69431, 0.744423, \
0.817768, 0.991047, 0.886702, 0.283607, 0.46127, 0.808971, 0.0370969, \
0.804104, Indeterminate, 0.631843, 0.989403, Indeterminate, 0.35993, \
Indeterminate, 0.960504, 0.749909}
*)

• Thank you for the response. I have to admit that I did not understand it. Could you please kindly explain, what is ElementMeshPlot3D  intended to do? The second part of your answer I unfortunately do not understand at all. Even not what did you try to achieve and what did you achieved as the result. – Alexei Boulbitch Jul 15 '15 at 8:15
• @AlexeiBoulbitch ElementMeshPlot3D plots the the surface with the heights given by the first argument over the element mesh in the second argument. See mathematica.stackexchange.com/questions/57366/…). I'll think on the second part -- got to go.... – Michael E2 Jul 15 '15 at 19:31
• @AlexeiBoulbitch Normally I use NDSolve in the form slv = NDSolve[..., u,...] so that slv is an InterpolatingFunction[...]. I find the functions easier to use. You had slv = NDSolve[..., u[x, y],...] so slv is InterpolatingFunction[...][x, y]. Thus Head[slv] is just InterpolatingFunction[...]. I could have rewritten the NDSolve code, but I thought this would be easier for folks to copy and paste. In any case Head[slv] evaluates the solution on whatever is given to it as arguments. – Michael E2 Jul 16 '15 at 0:18
• @AlexeiBoulbitch Yes, the second part makes me think the interpolation is buggy. But I wish I knew why reevaluating the same function on the same inputs gives different results. That makes me wonder if I'm doing something wrong. I too was wondering if the number of holes was a factor. -- Oops, I just realized I forgot to include the NDSolve with the Indeterminate extrapolation value. – Michael E2 Jul 16 '15 at 13:55
• @MichaelE2, no you are not doing anything wrong. It has to do with caching; it's complicated and I'll write more once I have a better understanding of it. But it's a bug. – user21 Jul 16 '15 at 14:17