Profiling Results + Proposed Optimizations

[dseyler]

Profiling Results

Over the last week, I profiled our solver and looked for places in the code that could be optimized for runtime efficiency. Overall, these optimizations gave a speedup of roughly 4x for my test case. Most of these savings came from optimizing assembly.

Rather than open ~10 different issues, I figured it might be more useful to give a brief summary of each optimization and discuss which improvements would have the biggest impact and least risk of complication. Also, I think it would be great to discuss how to make profiling a routine part of maintenance going forward.

Test case:

• cardiac mechanics 'struct' simulation coupled to svZeroD
• 584k linear tetrahedral elements
• Holzapfel-Ogden material model with prescribed active stress
• Trilinos MueLu

Profiling:

• Profiled on Apple M3 and Sherlock Skylake CPU (1 and 24 cores)
• Timed each phase (assembly, linear solver, etc.) with in-code timers (implemented in this branch)
• Static counters on Array, Array3, and Vector allocations
• Operation count / timing using perf @ 1999 Hz

Proposed Optimizations

Below is a brief summary of each optimization and the branch I implemented it in, ranked roughly in order of greatest speedup/least risk to least speedup/most risk.

Note that all of these features were implemented with heavy use of AI and definitely need to be checked over thoroughly and cleaned up before opening any PRs. That being said, they do pass all test cases and give identical Krylov counts and results for my profiled simulation.

• mat_fun optimizations: We already recently merged a PR using fixed-size Eigen arrays to template mat_mul. This can be expanded to other functions in mat_fun.cpp/h, including double_dot_product and dyadic_product. This may be better refactored into @zasexton's new Math module.

  • double dot product branch
  • dyadic product branch

• Templated viscous stress: Both potential and Newtonian viscosity models can be templated on nsd with fixed-size eigen maps being used within the functions, similar to how pk2cc already works.

  • templated viscous stress branch

• Struct3d stack buffers: struct_3d currently allocates 13 heap Arrays/ Vectors per call per Gauss point. As a result, allocation and destruction account for the majority of time spent in struct3d. This can be avoided by replacing these with stack-allocated arrays. The Array and Vector classes already have non-owning constructors, and one can be added for Array3 as well, although I think @michelebucelli pointed out this may be confusing for developers. Another possibility is to create a stack-allocated scratch buffer that is resized for each mesh and passed in as an argument to struct_3d(). This could also help with readability, as it would cut the number of arguments from 19 to 10. Maybe there's a better OO way to implement this change, so I'm open to suggestions.

  • struct3d branch

• struct3d linear element loop: Many of the quantities computed in struct3d depend on the derivative of the shape function, Nx and are computed once per Gauss point (4x for linear tets), but Nx is constant across the whole element for linear tetrahedral shape functions, so these Nx-derived quantities only need to be computed once per element. This gives a huge savings in Assembly for linear tetrahedral elements, but definitely needs to be treated carefully as I'm not sure how to best handle this if we allow multiple element types in one simulation. Currently, this is only implemented for struct3d but could be extended to other physics types too.

  • struct3d linear element loop

• bar_to_iso: bar_to_iso Is one of the most costly functions in struct simulations and can be sped up 2.4x using Eigen arrays. Note that this benefit is greater if the previously mentioned struct3d linear tet optimization is not implemented as it will need to run 4x more often. Claude also suggested another formulation: PP · CC _bar· PPᵀ = CC_bar − (1/n)·w·ciᵀ − (1/n)·ci·wᵀ + (cw/n²)·ci·ciᵀ where w = CC·c and cw = c·w. This version is 11x faster but is less readable than the more familiar PP : CC_bar : PP^T equation and requires CC_bar to be symmetric (I think this should always be the case, though @divya-adil ).

  • bar_to_iso branch
  • bar_to_iso with Claude's math

• Trilinos reuse: The Trilinos MueLu hierarchy is currently rebuilt every Newton iteration, and in my test case, this rebuild takes longer than the linear solve itself, negating the gains given by a better preconditioner. Ideally, the Trilinos hierarchy can be cached across time steps and only rebuilt when the tangent matrix has changed significantly such that a new preconditioner is needed. As discussed with @sujaldave, the safest option would be to rebuild it every time step rather than every Newton step, but other less strict criteria could be used, like a number of time steps or when the number of Krylov iterations starts to increase. It may even be possible that for some simulations, the Trilinos hierarchy never needs to be rebuilt, and in fact, Trilinos offers several different reuse levels that range from full reuse to complete rebuild.

  • Trilinos reuse

• All optimizations together: And finally, all optimizations on one branch

Let me know your thoughts on any of these ideas and which would be worth opening issues or PRs for! I think the first 3 + bar_to_iso would be fairly easy to merge and have little chance of collateral damage. The linear element loop is definitely the biggest win but would also need a lot of testing since an element that goes down the wrong path would not show any obvious problems.

[michelebucelli]

This is excellent @dseyler! Just to be sure I'm reading the graph correctly, the improvements are incremental, right? E.g., does the rightmost column correspond to all the optimization together?

[dseyler]

Yes! They're added sequentially. The improvements are somewhat dependent on the order in which the optimizations are applied