Add option to disable caching in FE::reinit()

[nmnobre]

Hi,

It turns out the problem I initially thought was to do with tetrahedral refinement is not specifically related to our refinement procedures, it is just more likely to happen when we refine.

There are cases in which two consecutive elements are sufficiently similar to trigger and reuse cached shape functions when, in fact, the orientation of the second element doesn't match that of the first. For example, two parallel edges aligned with the y-axis on neighbouring elements might not have the same orientation because the second has got a O(1e-16) difference in the x-coordinates of its two vertices for example. This is bound to happen when we refine and the resulting coordinates can't be exactly represented numerically (in base 2).

I've added a case to Vector FE Ex. 3: you expect grid_size=5 refine=2 to be the same as grid_size=20 (our default), but unless we disable caching, that isn't the case.

Of course an alternative to this would be, just like we do for caching, to do fuzzy comparisons for orientations as well. Opinions?

Cheers,
Nuno

[roystgnr]

Of course an alternative to this would be, just like we do for caching, to do fuzzy comparisons for orientations as well.

I don't see how this would work. There has to be some single point at which an orientation flips; it's not something that can be non-transitive like the similarity we use in caching. I'd love to find out how I'm wrong, though, particularly if it's a foundational sort of mistake on my part. At some point we need to come up with orientation rules that can work with mesh displacement, and I haven't thought of anything better than "use unique_id() when it's enabled".

My first thought was to start doing non-fuzzy exact comparisons for caching, but that won't work well: a cartesian grid which should cache cleanly would then start hitting false negatives with any rotation or even with most grid spacings.

My only other idea is that we could make the "fuzzy" tolerance controlled by the FE type. Using a nodal or discontinuous type, it specifies <= TOLERANCE*TOLERANCE, and you still get the maximum reasonable caching; using a type that depends on orientations for vector or high-order scalar shapes, it specifies <= 0, and you get safety.

[nmnobre]

There has to be some single point at which an orientation flips; it's not something that can be non-transitive like the similarity we use in caching.

So, similarity is a relation between two elements, and I agree it's non-transitive. The question that needs to be asked is: if two elements are similar, can we guarantee they have the same orientations. If two elements are similar, we know vectors (0) -> (m) and (0) -> (n), where m and n are vertices of some edge, in the two elements are similar, i.e. their difference has l1-norm <= tol. It's obviously the case this doesn't guarantee (m) -> (n), i.e. the edge, is similar by the same standard. But, by the triangle inequality, we know it's similar if we relax the constraint to <= 2*tol. So could we not define our lexicographic comparisons with this in mind? EDIT: this is probably irrelevant anyway, because we want to one day support displaced meshes (see below).

My first thought was to start doing non-fuzzy exact comparisons for caching, but that won't work well

Agreed.

At some point we need to come up with orientation rules that can work with mesh displacement

Since this is a worry, and I agree we'll need to use a unique node id, and since in that case our similarity test is not good enough, i.e. it's no guarantee the orientation is preserved, I'm tempted to go with:

My only other idea is that we could make the "fuzzy" tolerance controlled by the FE type. Using a nodal or discontinuous type, it specifies <= TOLERANCE*TOLERANCE, and you still get the maximum reasonable caching; using a type that depends on orientations for vector or high-order scalar shapes, it specifies <= 0, and you get safety.

If we decide to go with this, do we still keep --disable-caching or do I remove that?

EDIT: Or, we could instead define --enable-caching-for-orientation-dependent-elems (you get the idea) that the user could use to get extra performance once they know their mesh is amenable to caching...
Or, yet another option, for these orientation-dependent elements, we add an extra check on the caching, the second element needs to be similar AND preserve orientation.

[moosebuild]

Job Coverage, step Generate coverage on fa01248 wanted to post the following:

Coverage

	1fba56	#4119 fa0124
	Total	Total	+/-	New
Rate	63.32%	63.33%	+0.01%	100.00%
Hits	74333	74336	+3	10
Misses	43056	43048	-8	0

Diff coverage report

Full coverage report

This comment will be updated on new commits.

[roystgnr]

If we decide to go with this, do we still keep --disable-caching or do I remove that?

I definitely like having an option to disable caching. It's kind of like the option to use higher than double precision: we should never expect to need it, but if we ever suspect we need it then it's invaluable to be able to easily test that suspicion, either so we can quickly find and fix the problem or so we can rule out that suspicion. I guess we would want it to be a configure option rather than a command line option if we discovered that the extra runtime test hurt performance, but if we do the runtime option right then I'd be so surprised if someone saw a performance hit from it that I wouldn't trust their numbers without testing it myself. (I'll add some review comments re: doing it right)

the second element needs to be similar AND preserve orientation.

This isn't crazy. It'd be much more productive than the "set tolerance to 0 for orientation-dependent elements" idea.

[roystgnr]

This should probably be a static member - even if anyone is doing something as weird as turning caching on for some FE objects and off for others, that means they're deep enough in the weeds that they could figure out how to do that by manipulating the static member.

[nmnobre]

I've removed this member entirely. I can put it back if you think we should?

[roystgnr]

And if we make this a static member, could we initialize it with a singleton? I really hate doing any kind of string comparisons for variable identifiers in places where they're going to be hit more than once; I still remember the first time I profiled a hypersonics code that turned out to be spending a plurality of its time in string parsing...

[nmnobre]

I'm now (partially) using the Singleton functionality though I actually never create a Singleton object as it looks unnecessary here.