Anisotropic FEM forces float32 matrices

[VictorBarnes]

Question/Support Request

Hi LaPy developers,

I noticed the Solver class uses different default dtypes depending on whether the isotropic or anisotropic FEM method is used.

In Solver._fem_tria, the stiffness and mass matrices are each cast to a csc_matrix without an explicit dtype (taking the default float64). However, in Solver._fem_tria_aniso, the stiffness and mass matrices are each cast to a csc_matrix with dtype=np.float32. This means using the aniso method reduces the precision of the LBO matrices and the resulting eigenmodes and eigenvectors of the generalised eigenvalue problem also come out as float32.

This has become a problem because I have been performing large operations on the eigenvectors of a high resolution mesh and the float32 eigenvectors sometimes trigger matmul precision warnings. I'm also seeing NaNs appear in some of the outputs of these computations. However, these issues don't occur when using float64 eigenvectors.

Was there a reason for why the mass and stiffness matrices of the anisotropic LBO were cast to float32? If it is not essential to reduce the precision would you be able to remove the explicit dtype casting in _fem_tria_aniso? If it is essential, would you be able to add a dtype parameter to the Solver so the user can specify the precision they want?

Thanks in advance.

Environment

• LaPy Version: 1.5.0

[m-reuter]

Thanks for reporting. I don't think that there is a reason for explicit dtype setting. It is probably a leftover from some earlier version. Please feel free to create a pull request with the proposed changes. We can remove the explicit casting and make it similar to the fem_tria . Would be great if you can run a few tests with your use-case and, if possible, check other ways that this function is being called from within lapy to avoid that we have type incompatibilities anywhere in the code.

[m-reuter]

I actually created a PR for this (allowing users to specify dtype explicitly, default to float64 instead of getting it from the mesh). See #105 . This should solve your issue.

[VictorBarnes]

Thanks for the fix!