[timofeymukha]

Alrighty, this PR adds classes and routines to work with nmsh and re2 mesh files. The code is written by Fable 5.1, with a lot of steering from my side after its original implementation to try and fit the structure of the package better.

• the io module contains re2 and nmsh io, plus the corresponding Data classes. This corresponds heirarchically to pymech's hexdata.
• the datatypes module adds the corresponding types for the meshes, e.g. NmshMesh, both derived from a CornerMesh base type. The latter provides a converter routine to Mesh.
• The Router type is enhanced to work with records to shuffle around non-trivial numpy dtypes.
• A helper for linearly distributing data across ranks is added and used throughout the package, removing some duplication.
• Some tests.

This is the gist of it. I haven't verified every line of code, but I did go through each routine to see if it looks reasonable. Some routines in more details, like the one for building the trilinear basis (which is why it is so commented). Basically spent about a day on it, so not pure vibe coding :-).

[adperezm]

Nice! thank you.

I have just some comments that come from my own lack of knowledge perhaps.

With the functionality, I will unfortunately need to just trust the tests, as I have nowhere near enough time to check, given that I am not familiar with the formats in question (re2 and nmsh)

[adperezm]

This "record" terminology. Is it a standard?

I knew this as a "custom data type" and I believe that would be clearer to read. But ignore if this is standard.

I am bit confused here. Sorry if I am dense, so here you create one custom data type that has as many bytes as the full message so you just send one "count" (for lack of a better word) instead of >2**32 if you need to?

[timofeymukha]

Not sure if it really is a standard but record is kind of an alternative synonym to a struct, I'd say. Which is of course a bundle of binary data, which is what the record here is as well. The root may be in Pascal :P, at least there you write stuff like

	sphere = record
			radius: real;
			volume: real;
			surface: real
		end;

I remember learning it at university :D.

I think the idea is that you teach MPI to pass these around and then the limit of 2^32 applies to the structs. Since each packs quite a bit of data, the effective limit is then much larger.

[adperezm]

Yes but I still dont get it.

Here rec_t = MPI.BYTE.Create_contiguous(np.dtype(dtype).itemsize)

You are creating a new dtype that has an itemsize of the dtype that you pass. So in the end, if you pass an array of doubles, then you create a new dtype that has also size of double, so how is the count different?. Are you maybe using this somewhere where you pass a dtype that is different?

I have used this type of thing before, when I was implementing RMDA for interpolation:

pySEMTools/pysemtools/interpolation/interpolator.py

Lines 202 to 204 in e2cd9f4

	row_dtype = MPI.BYTE.Create_hindexed_block(blocklength=bytes_per_row, displacements=row_displacements)
	row_dtype.Commit()
	# Send

So what I would expect here is that you create a datatype where "one count" is more bytes than you original data type, such that you sort of reduce the total count below the int32 limit. But I do not see how that would happen here.

Again, I might just be too dumb haha.

[adperezm]

Sorry if I am annoying, I am just trying to understand a bit, because I do not want to have complex functions for the sake of complexity. If there truly is a need for something we keep it but I do not like overengineered solutions.

[adperezm]

Maybe to clarify, if I have an array

array = np.ndarray(2**33, dtype=np.double)

That has more items than int32.

If I just do (as is done in the code now):

        if key not in self._record_types:
            self._record_types[key] = record_datatype(array.dtype)
        rec_t = self._record_types[key]

That rect_t will just be a dtype of the same size of array.dtype

The way to overcome the int32 limit would be to create a new dtype that has more bytes than double, such that the count reduces, so one would need to do something like:

block_dtype = np.dtype([
    ("data", np.float64, (8,))
])

So we pack 8 doubles into one count for the new dtype

then rec_t = record_datatype(block_dtype) would solve the int32 problem. But right now I do not see anything like that.

Or perhaps I am wrong with this.

[timofeymukha]

I think you are right that if you pass an array of doubles nothing happens, and this is intentional because MPI already knows how to handle that, right? This thing is added for passing around mesh entities like

ZONE_DT = np.dtype(
    [
        ("e", "<i4"),
        ("f", "<i4"),
        ("p_e", "<i4"),
        ("p_f", "<i4"),
        ("g", "<i4", (4,)),
        ("t", "<i4"),
    ]
)

There are a some more defined in nmsh_file and re2_file. So one ZONE_DT element is the "record" and now MPI can pass these around.

[timofeymukha]

I could also be totally wrong though, and you worked with this stuff far more than I have! I basically thought this is like the stuff in mpi_types.f90 in Neko.

[adperezm]

Okay! I think I see now. array.dtype is used and it is some form of custom data type. We need a new mpi dtype that can communicate it.

I understand now. I think I got sidetracked by some comment I saw in the code regarding the int32 limit.

Using this would not fix the int32 limit by itself. The user must have already packed their data in a different way for that. But this change does allow to have custom data types, since before we were constrained to the basic ones.

Thanks!

[timofeymukha]

I think the comment may have appeared because I poked the LLM about what that was and it decided to focus on that in the comments :D