Generic implementation of a multicomponent Datatype class

[tlunet]

To follow up on the discussion started in https://github.com/Parallel-in-Time/pySDC/commit/04fdedc56e52258fa88ea98aae5602d5944fc7c0#r138624434, @brownbaerchen @lisawim (and also centralize all that have been done on this topic, don't hesitate to complement the description ...)

Workspace dev branch : https://github.com/brownbaerchen/pySDC/tree/multicompmeshv3

Current testing : https://github.com/brownbaerchen/pySDC/blob/multicompmeshv3/pySDC/tests/test_datatypes/test_multicomponent_mesh.py

Requirements :

• easy-to-use and generic implementation
• most operations (including numpy operators) conserve the type
• allow to access and modify the different components using a view on the same global array containing all components
• use for other current data classes : imex, multi-implicit, (particules ?)

Some preliminary works (some time ago ...) :

• base generic class : https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/playgrounds/datatypes/base.py
• application for mesh and others : https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/playgrounds/datatypes/implementations.py
• some tests for consistency : https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/playgrounds/datatypes/check_consistency.py
• some tests for performance : https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/playgrounds/datatypes/check_performance.py

$\Rightarrow$ main idea is to have a base multi-component class that can be used later for fenics, petsc, etc ...

[brownbaerchen]

To summarise the state of the development branch: We have an implementation that derives from Mesh and hence is only compatible with numpy stuff at the moment. What it does is expand the mesh along an added first axis for each desired component and returns a view onto the component when requested with .component. Making a specific multicomponent mesh is as simple as

class imex_mesh(MultiComponentMesh):
    components = ['impl', 'expl']

The implementation is a bit hacky and there are a few minor pitfalls. But with a bit of documentation, I believe it works well enough. In particular, with the above imex_mesh, all tests pass.

So we satisfy most of the requirements already. The main question is how far this should go. We can extend this to the particle datatype, but it's not as easy as for imex because of the mass that is stored along side the position and velocity. It's not difficult to do, but it does require overloading ufuncs and some care. Overall, I don't know if it would make the implementation much simpler. Also, I have not used things like fenics and petsc. Since this implementation is linked to numpy, I don't think it's much use there, but maybe someone else should comment on this.

So, to conclude: Should we merge the current state now so that @lisawim can use it for the DAE stuff and worry about how we can streamline some other stuff with this later? Or does anybody want this to be more thorough now?

[pancetta]

Can you also replace the comp2_mesh, too? I always wanted to merge all of these numpy-based data types into one core class, with whatever components the "users" want. So, I'm a big fan of this!

Not sure you need to include the particle data type, at least not directly. It's ok to have this around for a little longer.

[tlunet]

I would like to have a look at it during the week-end, especially since I'm not sure we still need to overload the ufuncs from numpy ... also, there may be a way to be more generic and avoid overloading of __getattr__ using metaclass (see for instance the preliminary work in the datatypes playground).

So you can merge it already if @lisawim wants to work with those classes for the DAE sweeper (could also allow to see other non anticipated pitfalls), and just let this issue open so we can continue exploring. Well defining a generic multi-component dataclass could allow to have generic implementation of the eval_f, solve_system and integrate for all the problem classes, see https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/playgrounds/datatypes/check_consistency.py.

End goal, this would greatly ease the usage and development of pySDC for additional problems (in particular for benchmarks ...)

[brownbaerchen]

Can you also replace the comp2_mesh, too? I always wanted to merge all of these numpy-based data types into one core class, with whatever components the "users" want. So, I'm a big fan of this!

That's already done on the development branch.

[pancetta]

In any case, if we do this and if this actually replaces the IMEX mesh, we'd probably need more documentation to the new classes.

[lisawim]

Does someone say "documentation"? At this moment I'd be fine with merging the DAEMesh, and I agree that we need a detailed documentation for the multi-component stuff here for more clarity about how it exactly works.

[brownbaerchen]

I am wondering where we should stop using this. So far, I was only thinking about integration-specific labels for the components, but why stop there? With this it is very easy to put tailored meshes in problems. For instance, in van der Pol, we have variables $u$ and $v=u_t$. Why not do the following?

class VdP_Mesh(MultiComponentMesh):
    components = ['u', 'v']

class vanderpol(ptype):
    dtype_u = VdP_Mesh
    dtype_f = VdP_Mesh

I don't really care for doing this now. But I think we should consider this when adding problems in the future. However, when having multiple implementations regarding integration, this is dangerous. I am not sure if we can "stack" the current MultiComponentMesh to do, for instance, an IMEX splitting on top. Anyways, feel free to ignore this rant.

[pancetta]

I like that! Yet another dirty code where I did not find a better way to handle this. This would be great to have!

[lisawim]

Such stacking of components could be very interesting for different cases, I like the idea!

For DAEs, it would also be interesting to stack the different components. Thinking of SDC for semi-explicit DAEs, the right-hand side can also be separated into implicit and explicit parts, therefore we would get something like

f.diff.impl, f.diff.expl, f.alg.impl

that requires stacking of diff and alg components. For the future this could be nice when developing a corresponding sweeper. However, this idea can wait for now.