#22 Field.evaluate: uses arbitrary cell for evaluation, within tolerance (1e-3 by default)

[Václav Šmilauer]

For an unknown reason, Example 05 sometimes fails (about 80%) and sometimes passes. The field value at t=0.1 is 727.55 (mostly, fails) and 728.13 (sometimes, passes). The behavior is the same for py2 and py3. Time values are the same for both cases. It is the only (ATM) outstanding issue for tests being consistently passing at http://travis-ci.org/mmp-project/mupif .

[Václav Šmilauer]

I am not able to reproduce this recently, having run Example05 100x with python2 and python3, all results were correct. But it sees that there was no code update recently which could have fixed this. Is someone (Vít?) still seeing the issue?

[Václav Šmilauer]

Which reminds me of how debugging works:

[Borek Patzak]

I would suggest to print some additional information. For example, the identified source element, its vertex values, etc. In case the test fails again, the additional info would be extremely helpful. I suspect, for example, that the point to be evaluated is very close to the edge between elements (or vertex), and the field value is "sometimes" obtained from the other element.

[Václav Šmilauer]

Just compared field values (through saving to HDF5) in successful and failed example, and both values and meshes (including vertex coordinates) are exactly the same (in terms of float comparison).

The problem is in Field._evaluate, where there is a loop over all cells containing the point (with tolerance, which is 1e-3 by default (!), which happens to equal the cell size in Example05) and the value of the first cell intersecting the tolerance box will be returned.

This might be incorrect even for interpolated vertex-based fields, like in the attached image: if cell1 comes first in the loop, cell1's interpolation function will be used to evaluate the file in the point which actually lies in cell2 (see attachment for illustration).

For cell-based fields, where there is no interpolation done, cell value which comes first is used as the result. In Example 05, those are all matching cells:

Value for cell w/ LL corner at (0.0, 0.001, 0.0):  (727.55,)
Value for cell w/ LL corner at (0.0, 0.0, 0.0):  (727.55,)
Value for cell w/ LL corner at (-0.001, 0.001, 0.0):  (727.55,)
Value for cell w/ LL corner at (-0.001, 1.7347e-18, 0.0):  (727.55,)
Value for cell w/ LL corner at (-0.001, 1.7347e-18, 0.0):  (728.13,)
Value for cell w/ LL corner at (-0.001, 0.001, 0.0):  (728.13,)
Value for cell w/ LL corner at (0.0, 0.0, 0.0):  (728.13,)
Value for cell w/ LL corner at (0.0, 0.001, 0.0):  (728.13,)

Then depending on which cell comes first when iterating over cells, the value is either 727.55 or 728.55.

So what I propose is to make the evaluation 2-pass: first use zero tolerance (whatever eps is); if that succeeds, return that . If no cells match (due to rounding and the point is exactly on the cell boundary or something like that), use eps.

It still can be a problem if eps is too large; generally it is good to set the defaults relative to the current scale of the problem, e.g. by making eps negative by default (-1e-6) signifying that it is relative to some problem-specific dimension. Field could, for instance, ask its mesh for the minimum element edge size (and remember it forever) which should work well as some characteristic dimension.

[Václav Šmilauer]

I fixed Example05 by explicitly setting tolerance to 0.0 and evaluating the field in 1e-5,1e-5,1e-5 instead of 0,0,0 (which is exactly on the vertex). The bug was retitled since the underlying issue is not settled, although Example05 should work now without issue.

[Borek Patzak]

Good.
If I remember right, the problem with this specific example is that the field is cell-based, so picking the different elements can make a big difference.
In a case of interpolated, vertex-based fields, the effect of picking different element should have rather small effect (as all elemens sharing the vertex should return the same value at the vertex), depending on the value of eps.
The question is now if we should opt for proposed two-stage evaluation authomatically, or left it on the user. I would prefer to have simple well defined functions, without the complex logic. So I wote for leaving this on the user.

[Václav Šmilauer]

One option is to require the eps parameter to be passed to evaluate (not providing any default value). Obviously, it is a nuisance, but better than giving some arbitrary value which might cause results which are incorrect.

We can also set eps=0.0 by default, which is the only non-arbitrary value for absolute distance tolerance. In that case we risk that evaluation will fail due to holes in floating-point inside tests in geometrically contiguous cells.

From these two, I see the first as the lesser evil.

Default non-zero relative tolerance would be the best, but field would need to have an idea about its cell size.

[Vit Smilauer]

I can confirm that using default value for eps=0.001 led sometimes to wrong result
Field value at position (0.0, 0.0, 0.0) is (727.55,)
2016-03-30 08:37:49 ERROR:Example05.py:49 Test FAILED

The last commit set up eps=0.0. It works for all 100 cases

Field value at position (-1e-05, -1e-05, -1e-05) is (728.13,)
2016-03-30 08:39:06 INFO:Example05.py:50 Test OK

[Václav Šmilauer]

I added automatic padding with the default value of relPad=1e-5. This will add extra 1e-5*bbox_size[i] to all 6 bbox directions (along each axis). Thus the padding is relative to cell size (scaling automatically with the problem) and non-zero. The commit is https://github.com/mmp-project/mupif/commit/1c74cd255bce025778237df7fdbc61834504ad61 . Along with that, the tolerance value was set to 0.0 by default, it should not be needed at all.

I am quite happy with this solution, since it avoids finite-precision issues we've seen before, puts no additional burden on the user and scales automatically. I am open to comments you might have.

[Ralph Altenfeld]

That's good news for me.
The files which are evaluated in example 5 are from our project test case. I circumvented the problem so far by evaluating at a z-coordinate 0.0 + eps.

[Borek Patzak]

I have changed the way, in which cell based field is evaluated. When a point belongs to multiple cells (when its on shared edge or even in a vertex) then the returned value is obtained as average over shared cell values.
The old way was to simly return the value of the first cell found countaining the point. In same cases this could cause the problem reborted above.
But overall, I believe, that the new approach is much better form numerical and physical point of view.