Track regions that fields are valid over

[dschwoerer]

Resolves #2295

I am not sure it is particular useful for 2D - but at least the patches applied cleanly.

Could be extended, to give a speedup in more cases, as with the region tracking, boundaries can be excluded if they are not needed, at the cost of more lookup-maps (that need to be computed once).

I implemented this only for Field3D because otherwise the 2D-3D look-ups are needed, and I only needed it for Field3D.

[ZedThree]

Please could you give a high-level overview of what this is intended to do, as well as how it works? I gather it's for making sure operations are done over consistent regions, but it looks like it does the operations over the intersection of regions. Do we want that behaviour over just checking the regions are consistent?

I think you also need to add docstrings on all the new functions, as well as some low-level comments. I found Mesh::getCommonRegion quite difficult to follow, for instance.

I'm also a bit worried about performance. There's now more map lookups for everything arithmetic operator, even after the initial creation of the lookups.

Please could you add some before/after timing comparisons, at least for an optimised build?

[ZedThree]

Suggested change

	if (this->size() != other.size()) {
	return false;
	}
	for (auto i1 = this->begin(), i2 = other.begin(); i1 != this->end(); ++i1, ++i2) {
	if (i1 != i2) {
	return false;
	}
	}
	return true;
	return std::equal(begin(), end(), other.begin());

Should probably be implemented as a free function rather than a member

[dschwoerer]

std::equal in the 3 argument version is actually missing the size check:
https://en.cppreference.com/w/cpp/algorithm/equal

[ZedThree]

Good spot, the four argument version should do the trick

[dschwoerer]

Please could you give a high-level overview of what this is intended to do, as well as how it works? I gather it's for making sure operations are done over consistent regions, but it looks like it does the operations over the intersection of regions. Do we want that behaviour over just checking the regions are consistent?

The idea is that each field can store over what region it has valid data. So if a field is initialised from a constant value, it is valid everywhere. A derivative would only be valid in the region without boundaries. If you multiply it with a derivative, it is still only valid in the region without boundaries, thus the intersection is used.

I think you also need to add docstrings on all the new functions, as well as some low-level comments. I found Mesh::getCommonRegion quite difficult to follow, for instance.

Sure, will do 👍

I'm also a bit worried about performance. There's now more map lookups for everything arithmetic operator, even after the initial creation of the lookups.

Thinking about it again, the map can be avoided, and instead a flattend "2D" array could be used.

Please could you add some before/after timing comparisons, at least for an optimised build?

I guess that would depend on the size of the meshes. It should be particular bad for a small mesh. I tried the example/performance/arithmetic for benchmarking, and am a bit confused about the result (all with --enable-optimize=fast and no guard cells)

with map on a mesh with 2 * 2 * 16 with "RGN_NOZ" and "RGN_ALL" set:

TIMING	minimum	mean	maximum
Fields:	0.536 us	0.572 us	11.163 us
C loop:	0.046 us	0.053 us	0.132 us
Templates:	0.201 us	0.209 us	0.681 us
Range For:	0.138 us	0.163 us	10.706 us

map-free on a mesh with 2 * 2 *16:

TIMING	minimum	mean	maximum
Fields:	0.689 us	0.721 us	11.445 us
C loop:	0.026 us	0.037 us	0.235 us
Templates:	0.184 us	0.200 us	0.976 us
Range For:	0.127 us	0.146 us	0.403 us

map-free on a mesh with 2 * 2 *16 with "RGN_NOZ" and "RGN_ALL" set:

TIMING	minimum	mean	maximum
Fields:	0.569 us	0.587 us	1.570 us
C loop:	0.047 us	0.053 us	0.178 us
Templates:	0.208 us	0.218 us	0.762 us
Range For:	0.141 us	0.159 us	0.401 us

current next on a mesh with 2 * 2 * 16:

TIMING	minimum	mean	maximum
Fields:	0.660 us	0.767 us	12.847 us
C loop:	0.027 us	0.040 us	0.160 us
Templates:	0.180 us	0.198 us	0.619 us
Range For:	0.130 us	0.144 us	0.562 us

[ZedThree]

Oh that arithmetic performance example is maybe a little out of date now. Timings for a full model might be more instructive, like blob2d or elm-pb.

[dschwoerer]

I tried the blob2d model, but the variation is to large to see any differences.

I am however thinking of extending this, and also setting the region after taking a derivative the region it is valid for. That also requires some code to extend the region, e.g. after communication or appling boundary conditions. However, that would be incompatible with some models, as BCs implemented in different ways would not be captured automatically.

That could allow to scale a bit better towards small grids, as the halo would not always be included in arithmetic operations, but that would require more work - so not sure it is worth it ...

[bendudson]

Thanks @dschwoerer ! I think this looks pretty good, but have a couple of questions:

1. Is this something that should always be done, or is it a correctness check that could be disabled if CHECK == 0 for example?
2. If a user performs an operation "by hand", e.g. outerloop calculations, what happens if they don't set the valid region?

Happy new year!

[dschwoerer]

Thanks @dschwoerer ! I think this looks pretty good, but have a couple of questions:

I am happy to also add this to docs, but I am not sure where.

1. Is this something that should always be done, or is it a correctness check that could be disabled if `CHECK == 0` for example?

I think it should also be done. Right not setRegion() is not called in a lot of places, but if it is, it can result in performance improvements, as e.g. boundary regions can be skipped in the computation.

2. If a user performs an operation "by hand", e.g. outerloop calculations, what happens if they don't set the valid region?

As long as setRegion() is not called, nothing changes.

It is not called in a lot of cases, as I was not sure how to handle users writing data. For example, the derivatives could call setRegion() with the region that they have computed the field, but then if the user implements some custom boundary function, the region will not be updated, and things might break.

Maybe in the future we can add this as optional feature, to call setRegion() for the derivatives, and for the build-in boundaries extend that region appropriately. For outerloop it should be fairly straight forward, as the user does everything, and can ensure that if setRegion() is used, it is used correctly.

Happy new year!

Happy new year and thanks for the questions :-)

[ZedThree]

Is there another way to write this without the #ifdef guards and repeated conditionals? Maybe BOUT_OMP(single)? If not, then I think this needs a good comment to ensure it doesn't get accidentally removed

[dschwoerer]

I added comments. Please let me know if there is still anything unclear.

[bendudson]

Bikeshedding perhaps, but I think the name getDefaultRegion could be confusing: It's getting the field region, using the given region name as the default, rather than getting the default region. Perhaps getRegionWithDefault?

[bendudson]

Since getRegion gets the region with the given name or ID, perhaps getValidRegionWithDefault is more explicit that it's the valid region that's being requested. Too long?

[dschwoerer]

I changed it to getValidRegionWithDefault - that is probably easier to understand. The length seems ok to me 👍

[ZedThree]

This explanation could do with expanding. I think I now understand what's going on, but this diagram is definitely confusing by itself

[dschwoerer]

I tried to explain it, but apparently failed to do so.
Would you mind adding a comment, as you understand it now? If not, I can try again, but I am not quite sure how to do it better :-(

[ZedThree]

This function finds the ID of the region corresponding to the intersection of two regions, and caches the result. The cache is a vector, indexed by some function of the two input IDs. Because the intersection of two regions doesn't depend on the order, and the intersection of a region with itself is the identity operation, we can order the IDs numerically and use a generalised triangle number: $[n (n - 1) / 2] + m$ to construct the cache index. This diagram shows the result for the first few numbers.

These indices might be sparse, but presumably we don't expect to store very many intersections so this shouldn't give much overhead.

After calculating the cache index, we look it up in the cache (possibly reallocating to ensure it's large enough). If the index is in the cache, we can just return it as-is, otherwise we need to do a bit more work.

First, we need to fully compute the intersection of the two regions. We then check if this corresponds to an existing region. If so, we cache the ID of that region and return it. Otherwise, we need to store this new region in region3D -- the index in this vector is the ID we need to cache and return here.

[github-actions]

clang-tidy made some suggestions

[dschwoerer]

Failure seems to be related to the fact that the --cflags do not include -std=c++17

I am not sure to what extend that was already an issue in the past while compiling against bout++, but it certainly is now required if you include mesh.hxx or field3d.hxx

[dschwoerer]

36f152d might be worth backporting to master?

[ZedThree]

We might need to revisit this as this is compiler-dependent

[dschwoerer]

👍 I am not sure how widely that is used, thus I would wait whether anybody ever complains ...

[ZedThree]

Suggested change

//#include "bout/output.hxx"

[ZedThree]

This function finds the ID of the region corresponding to the intersection of two regions, and caches the result. The cache is a vector, indexed by some function of the two input IDs. Because the intersection of two regions doesn't depend on the order, and the intersection of a region with itself is the identity operation, we can order the IDs numerically and use a generalised triangle number: $[n (n - 1) / 2] + m$ to construct the cache index. This diagram shows the result for the first few numbers.

These indices might be sparse, but presumably we don't expect to store very many intersections so this shouldn't give much overhead.

After calculating the cache index, we look it up in the cache (possibly reallocating to ensure it's large enough). If the index is in the cache, we can just return it as-is, otherwise we need to do a bit more work.

First, we need to fully compute the intersection of the two regions. We then check if this corresponds to an existing region. If so, we cache the ID of that region and return it. Otherwise, we need to store this new region in region3D -- the index in this vector is the ID we need to cache and return here.