Use geometry type independently#4275
Conversation
|
It is more subtle than described. The UFCx interface supports a single real scalar type. There is no point mixing, for example, in a straightforward way a double finite element and float geometry because the the float get promoted to double in multiplication, etc. It wasn't half-baked; it was deliberate (but there are of course other ways too).
|
Right. But why should we limit ourselves to what UFCx/FFCx support at this moment? We've been allowing custom kernels for years. This is still a draft and I am testing downstream for a combination of scalar_t = float16, geom_t = float64. Think of: I want dolfinx to store a precise float64 mesh, do the assembly loop, but I'll supply my own kernel that does the magic and handles the mixing of geometry of FE fields. And we can improve the UFCx/FFCx do in these cases in near future. |
There was a problem hiding this comment.
This looks like nice work towards #3543!
I do not like that this is a global registry: any change in a type needs to be registered not only in the export/wrapper, but also here. Maybe move the lookup table as static member into the python wrapper class? There any newly exported type needs to be made available anyhow. So for Function for example add static Function.cpp_registry, then a general (centralised, but stateless) cpp_type(python_type, dtype) can be implemented that does the lookup (based on the static registry attached to python_type).
There was a problem hiding this comment.
Good point, I've moved the registry to be local to each class.
There was a problem hiding this comment.
Is the scalar- and geometry type tuple the right general abstraction? Our exported types currently do not differentiate this, say for Function we have the export Function_float64. Also some types will never depend on such a concept, say for example Vector. I had in mind to just do a dtype based lookup. In case this is not flexible enough (in the future), specialised implementations for these types could be provided.
There was a problem hiding this comment.
I don't fully understand what do you mean. No, this is not a general abstraction, and see my comment #4275 (comment) Some classes will need more, some less types.
With this PR this is possible, but the naming convention is up to you,
dolfinx_wrappers::declare_objects<float, double>(m, "float32_float64");
dolfinx_wrappers::declare_la_objects<std::float16_t>(m, "float16");There was a problem hiding this comment.
Then the interface (key type) to the registry is not what we want. How about tuple[dtype, ...]? Any export is associated with a variadic length list of types associated to the template types in c++. That we can use as a general interface for type look up.
There was a problem hiding this comment.
I've moved registries and respective registry functions to be all local to classes, classmethods. So they have the right signature now, based on how many dtypes they require.
|
Updated now, ready for review I think. This is a minimal change and it allows independent use of scalar type and geometry types. I tested e.g. scalar type fp16 and geometry type fp64 downstream. This is not the whole story, as @garth-wells says. In the next PR I'd make also the kernel output tensor type independent, std::function<void(V*, const T*, const T*, const U*, const int*, const std::uint8_t*, void*)>which is the more desired approach, e.g. inputs in low precision and accumulate into higher precision. |
| (np.complex64, np.float32, "complex64"), | ||
| (np.complex128, np.float64, "complex128"), | ||
| ): | ||
| register_cpp_type(Form, _scalar, _geom, getattr(_cpp.fem, f"Form_{_name}")) |
There was a problem hiding this comment.
Just write those out in the class member dict directly? The typing hints of _cpp_object anyhow list them already.
schnellerhase
left a comment
There was a problem hiding this comment.
We should add a test case for these new templates. Could add a cpp unit test or extend the custom_assembler demo.
I don't plan to do it in this PR. Dolfinx doesn't do anything that would require testing these mixed precision types. And matched precisions are exported and tested. Once this is more explored in derived libraries, we can move some exports incl. their tests to dolfinx. |
… (HTTP 429)" This reverts commit afb3038.
Current use of geometry type different from Function type is half-baked. This PR uses independent type for geometry consistently.
In some places, there is the default behavior to deduce the geometry type as
scalar_value_t<T>from the Function typeT.Remember there are two concepts:
dolfinx::scalar T, for Function-like data, could bestd::complex<std::float64_t>, andstd::floating_point Ufor geometry data, cannot be a complex type.AI assistance: I used Claude Code Opus 4.8 to draft parts of this PR. I reviewed, edited, tested, and take responsibility for the final contribution.