Skip to content

Commit

Permalink
Preparation for PDAE stuff (?)
Browse files Browse the repository at this point in the history
  • Loading branch information
@lisawim
lisawim committed Feb 14, 2024
1 parent be4f248 commit 04fdedc
Showing 1 changed file with 10 additions and 2 deletions.
12 changes: 10 additions & 2 deletions pySDC/implementations/datatype_classes/MultiComponentMesh.py
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
import numpy as np

from pySDC.implementations.datatype_classes.mesh import mesh


Expand All @@ -13,9 +15,15 @@ class MultiComponentMesh(mesh):

components = []

def __new__(cls, init, *args, **kwargs):
if isinstance(init, tuple):
def __new__(cls, init, val=0.0, offset=0, buffer=None, strides=None, order=None, *args, **kwargs):

This comment has been minimized.

Sign in to view

Copy link
@brownbaerchen

brownbaerchen Feb 14, 2024

Contributor

why do you need this? The previous implementation was doing the same but more generic because it inherits this stuff from mesh. I just meant testing different shapes in init[0] in the test. But the implementation of the MultiComponentMesh should be able to cope with both PDEs and ODEs if you do sth like here

This comment has been minimized.

Sign in to view

Copy link
@lisawim

lisawim Feb 15, 2024

Author Collaborator

What is self.flags['OWNDATA'] used for? I wrote a test for adding two meshes mesh, mesh2 via mesh += mesh2 and attributes are lost again after addition.

This comment has been minimized.

Sign in to view

Copy link
@brownbaerchen

brownbaerchen Feb 15, 2024

Contributor

Hm... this flag is set by numpy and tells you if you're dealing with a view onto some object which owns the memory or if the object owns the memory itself.
For instance, mesh.alg is a view and does not own the memory. We don't want to add .alg and .diff to mesh.alg (i.e. get mesh.alg.alg) when doing ufuncs on this. That's why added this. But I was afraid that this would yield some side effects I did not consider...

This comment has been minimized.

Sign in to view

Copy link
@lisawim

lisawim Feb 15, 2024

Author Collaborator

Ok, I see and understand your intention! But I'd assume that only mesh passes the __array_ufunc__ function and not mesh.diff and mesh.alg as components. And mesh is then the view which always should have the components as attributes. So, does it make sense to remove the statement requesting if the memory is own or not and only use if type(self) == type(results)?

This comment has been minimized.

Sign in to view

Copy link
@brownbaerchen

brownbaerchen Feb 15, 2024

Contributor

Indeed, we cannot use this flag. However, I am not sure if the code works if we just remove it. Because doing mesh.alg + mesh2.alg with an ODE with only one component will try to assign mesh.alg.alg = (mesh.alg + mesh2.alg)[1] which does not exist.

I am currently trying this, which is basically doing the same, but only when accessing the components will it try to index things, thus circumventing the aforementioned issue. What do you think about this, @tlunet? It does feel a bit hacky. But all tests pass when implementing imex_mesh and comp2mesh derived from this, see here.

This comment has been minimized.

Sign in to view

Copy link
@tlunet

tlunet Feb 15, 2024

Member

First, you must call the parent __getattr__ and not raise an error in the last line of the overloaded __getattr__ method. Because currently, you simply remove access from all attributes specific from ndarray or mesh.

Second, I'm not so sure what's the problem here ... do we have the tests somewhere that describe what we want to be able to do with this new datatype ? (with the extended list of all operations ...)

This comment has been minimized.

Sign in to view

Copy link
@brownbaerchen

brownbaerchen Feb 15, 2024

Contributor

No. __getattr__ is a very funny method actually. It is only called when __getattribute__ "fails". So really it is called instead of raising an attribute error, but not every time you access anything.
In fact, the issue I see with this is the other way around. If you call a component something that is already defined by numpy, you cannot access it.

This comment has been minimized.

Sign in to view

Copy link
@tlunet

tlunet Feb 15, 2024 •
edited
Loading

Member

Oh yes indeed, I remember now ... guess you implemented it well then 😅. And to solve the mentioned issue, I guess a simple test in __new__ checking if the object already have an attribute of the same name (and raises an error if yes) would solve it.

This comment has been minimized.

Sign in to view

Copy link
@brownbaerchen

brownbaerchen Feb 15, 2024 •
edited
Loading

Contributor

No, I don't think that would do it. mesh.alg is just a MultiComponentMesh itself with no knowledge of its "parent", the big MultiComponentMesh including all the components. How would it determine if it needs alg itself or not? The check if its shape allows it to have .alg is not perfect. But at least it doesn't cause headaches if you use this class as intended. I guess you could include that in __new__. But I prefer the __getattr__ because otherwise you need it in both __new__ and __array_ufunc__ and maybe in __array_finalize__ as well, although I think not.

This comment has been minimized.

Sign in to view

Copy link
@brownbaerchen

brownbaerchen Feb 15, 2024

Contributor

Oh, nvd I understand now what you are saying! I was considering this. I just thought that I want to do as little as possible in these methods as they are called many many many times. So I'd rather write documentation telling people not to do this than enforcing it in the code every time someone initialises a mesh. But I know that there are benefits to enforcing this as well... So I don't know which is best.

This comment has been minimized.

Sign in to view

Copy link
@tlunet

tlunet Feb 15, 2024

Member

Can we have this developed in some dev branch in the main pySDC repo ? That way we could work together on it and find the best solution ... because currently, it's hard for me to see the exactly what are the problems and potential solutions.

Also, @lisawim do you have an issue somewhere describing all the elementary operations this datatype should be able to do (while keeping its original type and not becoming a numpy array again ...) ?

This comment has been minimized.

Sign in to view

Copy link
@brownbaerchen

brownbaerchen Feb 15, 2024

Contributor

I isolated the MultiComponentMesh from the DAE stuff in this branch. Feel free to push here! There are currently no unit tests for this in this branch, as Lisa has started specific to the DAEMesh datatype. But we can easily adapt those tests to this. Also, since all tests pass, I have good confidence that this works. Otherwise, this would mean IMEX stuff is really poorly tested in pySDC...

My suggestion is that we fix the generic datatype on this branch. Everybody loves solving merge conflicts. No problem. Then, we merge this into the main repo and then @lisawim can easily add the DAEMesh to the PR by merging upstream.
Are you fine with this? If @lisawim is in a hurry to merge her stuff, we can do that first and refine the datatype later.

This comment has been minimized.

Sign in to view

Copy link
@lisawim

lisawim Feb 15, 2024

Author Collaborator

With the current version of the MultiComponentMesh I have no issues there. Everything seems to work as it should be, i.e., all data type are the data types that are supposed to be. Only during writing a test I recognized the issue mentioned above.

I already started with writing a test, but we also can move this to later. This does not feel comfortable to me, but anyway. So, fixing the issue feels better to me, but here there is no hurry with that. This PR is already open for a while and can remain open even longer (at least for me). I'm therefore fine with your suggestion @brownbaerchen

This comment has been minimized.

Sign in to view

Copy link
@brownbaerchen

brownbaerchen Feb 15, 2024

Contributor

I added a simple test here. Do find anything is missing? I also checked that IMEX stuff is not slower or faster with this implementation.

This comment has been minimized.

Sign in to view

Copy link
@tlunet

tlunet Feb 15, 2024

Member

Maybe some assertion that the type is conserved after some operation like +=, **, np.exp, np.sin, etc ..., and the others that you already tested.

This comment has been minimized.

Sign in to view

Copy link
@lisawim

lisawim Feb 15, 2024

Author Collaborator

Looks good! Might it also make sense to test if the attributes are contained, i.e., testing

for comp in mesh.components:
    assert comp in dir(mesh), f"ERROR: Attribute {comp} not contained in the mesh!"

after applying +=?

This comment has been minimized.

Sign in to view

Copy link
@brownbaerchen

brownbaerchen Feb 15, 2024

Contributor

Looks good! Might it also make sense to test if the attributes are contained, i.e., testing

for comp in mesh.components:
    assert comp in dir(mesh), f"ERROR: Attribute {comp} not contained in the mesh!"

after applying +=?

I test this implicitly by calling asserts on the components after the operation. Note that the components are not in dir because they are not really attributes. I added a test that we can still access the components after calling np.exp on it.
if isinstance(init, tuple) and isinstance(init[0], int):
obj = super().__new__(cls, ((len(cls.components), init[0]), *init[1:]), *args, **kwargs)
elif isinstance(init, tuple) and isinstance(init[0], tuple):
obj = np.ndarray.__new__(
cls, (len(cls.components), *init[0]), dtype=init[2], buffer=buffer, offset=offset, strides=strides, order=order
)
obj.fill(val)
obj._comm = init[1]
else:
obj = super().__new__(cls, init, *args, **kwargs)

Expand Down

0 comments on commit 04fdedc

Please sign in to comment.