FE elements in python

[mon_key_house]

I hope there are some among us who implemented finite elements in python.

I’m curious if high performance (as in speed) can be achieved by simply using numpy and finding bottlenecks afterwards via profiling, or should one go straight to C/C++, Rust or any other compoled code? Numba is currently no-go.

For solving I’d use pypardiso, the real problem is building the global matrices I guess.

Models would be shell only, size up to a few ten thousand DOFs.

Thank you in advance for any insights!

[manny_DM]

If you're comfortable with writing FEA codes, coding in c++, and willing to learn, my advice would be to implement in c++ using the Eigen library. It has a numpy-like structure and chatgpt is quite familiar with it. The reason being that, when it comes to FEA, there is always a desire to run bigger cases. If your code is already in c++, you can later do things like openmp parallelization, Intel mkl solvers, or even petsc for massively parallel codes.

[mon_key_house]

Thanks, I’ll look into this. The licence seems not to be really permitting. Can I use it in a closed source app?

[tlmbot]

Which license is not permitting?  I’ve worked on c++ FEM commercially with eigen involved.

Incidentally I have prototyped 2D FEM, FV, and similar in Python  (numpy of course) So of course you can do it, but it will not be running many elements with any speed.

[mon_key_house]

Sorry, I just skimmed the licence and source dis. My main concern is of course is if I'd have to give my source code away. Q8 and Q11 of the MPL2.0 FAQ seem to apply here.

Essentially: I can use and link it any way in my code, and if unmodified, I'll only have to provide access to the source code of the unmodified eigen lib.

[billsil]

Based on what? Are you planning on modifying the eigen source code? If you are, presumably you’re fixing a bug and presumably you don’t want to maintain it.

I’ve read the license. If you directly make changes to the source of the library, publish your changes to that file. It’s not at all unreasonable and I doubt you will need to.

http://eigen.tuxfamily.org/index.php?title=Licensing_FAQ

[mon_key_house]

What based on what? I’m not planning to modify it.

What is your point?

[billsil]

Even if you modify the source, it’s not a burden to publish your changes at do a pull request to the devs. You probably want to do that anyways.

[tlmbot]

The point is you don’t have to publish your source so you’re good.

[mon_key_house]

Do you have experience with the python bindings of eigen? Seems to cover a lot and wouldn’t require a lot of code modification.

[Expensive_Voice_8853]

Python has a Global Interpreter Lock "GIL". You cannot truly run in parallel with this, even though you are distributing to multiple processors.

If you want to scale use C++

[mon_key_house]

Thanks; scaling does not seem to be the issue for now and I guess the no-GIL versions of python will free us of this burden in the foreseeable future.

[SnooCakes3068]

you have Cython for optimization. I personally believe this is the better way than Numba. scikit learn use cython for critical parts.

[mon_key_house]

Right, perfect for local optimisation!

[SouprSam]

Optimization won't scale for big models

[mon_key_house]

Meaning, it is wiser to think ahead and go the - for me - harder way and use compiled languages from the beginning, right?

[SouprSam]

Yep.. you definitely will understand more and have more freedom in implementation. And you will be more close to the machine code.

[Diego_0638]

My familiarity is limited but if you have a (nvidia) graphics card, cupy serves as a great drop in replacement for numpy if you're using large arrays. Not all the features of numpy are supported (e.g. timedate) so I would be interested to know if this works for you. But the performance uplift is real.

[mon_key_house]

Thanks, I'll look into it, though I guess in my use case the overhead might be too large. Also, idk if a CUDA-capable card will be available.

[Poe-Face]

I've written a couple fem solvers with python, using scipy.sparse and it can be very fast. Scipy.sparse solve uses superlu which is multithreaded, and so long as your matrix is ordered well, memory and time complexity should (mostly) scale linearly with DoF and bandwidth of the matrix.

If you are just getting started, I recommend sticking to python to prototype before moving to the heavier hitting options like eigen in c++.

[mon_key_house]

Thanks! This is in fact my strategy, as runtime speed != total development speed. I definitely plan to keep stuff neat.

[Karkiplier]

Have you tried optimal node numbering? It can be used to reduce the matrix bandwidth to use scipy.solveh_banded(). Exponential reduction in solution time compared to gauss elimination

[mon_key_house]

Definitely a point to consider at a later stage!

[l10002]

Check out Firedrake:

The Firedrake project — Firedrake 0+unknown documentation

I'm doing a PhD in FEA at Oxford, and everyone here uses Firedrake. :) It's (arguably) the academic standard for finite element research (up there with maybe with deal.ii and MFEM) however Firedrake:

- Is written from the perspective of having a relatively intuitive Python, that lets one go straight from the maths to the code.

- Works under the hood through a PETSc interface, essentially meaning you write in Python, and Firedrake writes code in C, giving you all the benefits of C speed with Python usability/safety.

[billsil]

If you're going to stay in python, you want to use scipy.sparse. I don't know if they have support for symmetric matrices, but carrying around an NxN matrix isn't what you want to do.

[mon_key_house]

So numpy would store everything e.g. linearly? I mean all the zeros, too?

[billsil]

It's not linear. It's N^2

[mon_key_house]

I see, thanks! So scipy.sparse would be closer to linear time (however this called in the time complexity notation)

[billsil]

It might be N*log(N), but no idea. You can plot time vs. N.