TL;DR : Writing an implicit Material Point Method solver. Would it be that much slower in Julia than C++ for a similar time-effort ?
Hello,
I am PhD student working on simulating creeping snow. I am currently planning on writing my own computational mechanics solver : a quasi-static implicit Material Point Method. For this purpose, I am wondering if you just better do it in C++ or in julia for speed reasons. Let me explain the structure and stakes of my solver.
The Material Point Method is a hybrid Eulerian-Lagrangian Method relying on a set of particles which represent the material, and a background mesh. It goes basically like this :
All the material informations and properties is stored on the particles
These infos are interpolated on the back ground mesh. This step implies a small loop (a dozen step to reach the surrouding grid nodes) within a bigger loop (over all the particles).
The equation of momentum is solved on the grid in the fashion of a Finite Element Method. Warning : since I want to simulate an implicit scheme, solving the mechanical equations involve inverting a matrix.
Info is interpolated back to the particles (similar way as step 2)
Particles positions and speeds are updated, grid is reset, ready for next timestep calculation.
So basically, the optimization challenges (to my understanding and knowledge) are the following :
managing big data structures (100k to a few million particles)
parallelizing as much as I can
managing big (parallelized loops)
operations on matrixes, including inversion
making my multiple dispatch efficient (to switch easily between small variants, 2D/3D, etc)
others ?
Here is my issue :
I am a new user to Julia and I love it, the readability, the multiple dispatch, the expressivity, the debugging (vs c++). Which would make my code much faster to implement than in c++ (which I am also rather new to).
However, this project is not just personal but professional. I really need/want to implement the fastest code I can, even if it takes more time. I would accept being up to 20% MAX slower in Julia than what would be reasonably achievable in c++, but cannot accept a 2 times difference (4 hours long vs 8 hours long simulation is a big difference).
Now, Iâve seen many threads saying Julia is 2 times slower, 30% slower, faster sometimes. And I understand it really depends on the situation. From what I understand, while Julia is slower in general, it would be easier for me to optimize the code in âsimpleâ Julia than in âcomplexâ c++.
So for an equivalent time-effort, would you think Juliaâs efficiency would be close enough to that of c++ for my purpose ?
I think the Julia code will probably be faster, because it is easier to profile, measure, and focus the attention of the performance on the bottlenecks (thatâs my experience vs. Fortran).
But even better is to take advantage of the available frameworks.
It really depends on your lvl on c++, you can do crazy performance critical decision in both language. However it will be easier to change the code to perf better in julia, c++ would require you to rewrite big parts to get a little perf sometimes.
If you choose julia I can only tell you to follow closely Performance Tips ¡ The Julia Language, make your code the best possible (time wize and allocation wize) before even considering parallel code, and only then you can go for it.
Doing it from scratch can be a good experience, however, donât think it will add anything to your phd (unless you wanna go back to the engenering path after it), so go for the tools out there and try to understand them.
If youâre not an expert in C++, I think itâs very likely that with 2 weeks of effort your Julia code will be faster than your C++ code. At my company we tried reimplementing a Python pipeline in several different ways, and within two weeks the Julia version was much faster than the others.
Basically, Julia can do most of the optimization-related things C++ can. You can parallelize pretty easily in fine-grained detail with Tasks, code is compiled via LLVM, itâs quite possible to write non-allocating code, you can easily examine generated LLVM or assembly, etc.
So the big difference is that with Julia you will get a working prototype much sooner, and be able to profile it in detail more easily. (You can of course do detailed profiling in C++, but itâs harder to set up and learn, compared to macros like @timeit in Julia.) Youâll be able to find the slow parts of your code sooner, and spend time optimizing them.
Thank you for your input, I will indeed keep the performance tips on an open tab at all times .
Contrary to what you said, I think writing my own codebase is extremely relevant for my PhD since it forces me to go back to the equations and derivation of the solver. And writing the code makes that i fully understand it and can modify it easily. All this makes that, in the end, I will have a deeper knowledge of MPM, mechanical models and computational mechanics in general.
Yes deriving the numerical method, do some numerical analysis and such is beneficial for your PhD, however while writing it from scratch is beneficial for you knowledge, Iâm not sure it is for the actual quality of the PhD thatâs all, I just donât want you to pass all your PhD coding which happens a lot and isnât really relevant. If your PhD is on applied math obviously. Coding should be like 60% o think in applied math world. Unless you have an industry behind the PhD wanting a new code.
Iâd choose Julia for the reasons others have mentioned above. However, the downside of Julia for me has mostly been that while writing working code is fast and itâs pretty easy to come up with elegant code, a lot of time is spent tracking down obscure allocations in performance-critical functions. The issues often come from using high-level language features such as closures etc. I probably shouldâve internalized the performance tips earlier to avoid too much headache. However, sometimes the allocations can pop up in pretty mysterious ways, or at least ways that are not easy to reason about. So itâs good to learn all the available tools like Cthulhu.jl and perhaps the new AllocChecks.jl.
Many performance related matters are not specific to Julia though, including issues related to parallelization or memory access.
Regarding using existing toolkits vs. writing your own, one approach is to write a functional implementation by yourself but once your code has reached some level of maturity, start looking into the established toolkits instead going down the path of adding features and improvements to your code. That way, you learn what writing a solver entails but donât end up spending thousands of hours on a project thatâll end up lying abandoned in Github once you move to other things. Even a minor contribution to an existing project is likely have a larger real world impact than trying to maintain your own.
I will definitely consider writing my solver in Julia, according to you, I will reach similar performance as c++ since optimizing will be easier, while my work will be more comfortable in Julia.
I did not know about MaterialPointSolver.jl. Itâs a great bonus, I would be able to follow an example in Julia with its optimization features.
It is still not clear why you would want to write a new solver, instead of using (or contributing to, if necessary) the existing one.
Donât get me wrong, it will probably be a great learning experience, but what you will be learning is programming which may not be closely related to your PhD topic.
If you topic is an application of this method, then writing a new package if there is an existing one is a rabbit hole and a waste of time you could be spending on research. If you topic is a new refinement of this method, then you may be able to modify existing software. In any case, you may want to talk to your advisor first about priorities.
The reason I want to implement my own solver is that MaterialPointSolver.jl uses an explicit scheme, where I want to use an implicit one. Let me explainâŚ
If we call Q and quantity involved in the momentum equation (speed, density, etcâŚ) and ^n and ^{n+1} designate two successive time steps.
An explicit scheme calculates the state at ^{n+1} only with the state at time ^n : Q^{n+1} = f(Q^n). This method is the most straightforward solution, benefits from fast calculations and is often more precise. However, it is only conditionally stable, the condition being on the maximum size of the time step.
On the other hand, an implicit scheme calculates the state at time ^{n+1} by solving a system involving the states at times ^nand^{n+1} : Q^{n+1} = g(Q^n, Q^{n+1}). Solving this system involves bigger time and allocation since we need to inverse a matrix. But this method is unconditionally stable, which means we can choose bigger time steps.
Using an implicit scheme is more convenient when dealing with long simulations and rather static or âslowâ deformations.
To the best of my knowledge, there does not exist an implicit, quasistatic, Material Point Method scheme for geomechanical purposes, with all recent efficiency and precision improvements of MPM, written in a fast, parallelizable language (Julia, C++, âŚ).
There lies the innovation of what I do, which will be valuable in my PhD. Plus, I also might need to create a new viscoplastic model, which will be easier to implement if I know everything about my code.
Thanks for the explanation. Given that the package is actively developed and appears to contain a lot of code besides the solver, I would still try approaching the maintainer about making a contribution first.
Donât get me wrong, nothing prevents you from developing your own package. Itâs just that you may end up duplicating a lot of effort for no good reason.
.