Templated viscous stress functions + linear shape function fix

[dseyler]

Current situation

Address issues #633 and #634, which involve templating viscosity models on nsd and skipping unnecessary computation when shape function gradients are constant across an element. These changes reduce assembly runtime by ~30% and fix a bug where shape function gradients were assumed to be uniform within wedge elements.

Release Notes

• compute_visc_stress_potential are now templated on nsd in mat_models.cpp
• Eigen views are implemented throughout viscosity models to keep intermediate array computations on the stack and condense several for-loops into one-line expressions.
• compute_visc_stress_and_tangent dispatches to templated models by nsd
• mat_symm and mat_dev are also nsd templated in mat_fun.h
• Svis, Kvis_u, and Kvis_v are allocated outside of the Gauss point loop, written over for each Gauss point, and only zeroed when no viscosity model is defined for an element (sv_struct.cpp, ustruct.cpp, fsi.cpp) or when viscosity arrays can be reused.
• recompute_visc flag is added to indicate when viscosity computation can be skipped (shape function gradient are uniform for all Gauss points within an element, so Kvis_u, Kvis_v, and Svis values can be reused).
• compute_visc_stress_and_tangent skips viscosity computation when recompute_visc == False
• default case added to compute_visc_stress_and_tangent so viscosity arrays are zeroed and do not carry the previous element's values when an element doesn't have a viscosity model defined
• struct_2d, struct_3d, ustruct_2d_m, ustruct_3d_m now require four new arguments: Svis, Kvis_u, Kvis_v, and recompute_visc
• compute_visc_stress_and_tangent now takes recompute_visc
• Fixed bug where lShpF should be False for wedge elements.

Testing

All test cases pass. Eigen operations were verified in isolation to be within machine precision of the operations they replaced.

Code of Conduct & Contributing Guidelines

• I agree to follow the Code of Conduct and Contributing Guidelines.

[claude]

Claude Code Review

This pull request is from a fork — automated review is disabled. A repository maintainer can comment @claude review to run a one-time review.

[codecov]

Codecov Report

❌ Patch coverage is 62.70492% with 182 lines in your changes missing coverage. Please review.
✅ Project coverage is 72.88%. Comparing base (c9838bd) to head (81edd1a).

Files with missing lines	Patch %	Lines
Code/Source/solver/ustruct.cpp	58.18%	97 Missing ⚠️
Code/Source/solver/sv_struct.cpp	56.89%	50 Missing ⚠️
Code/Source/solver/mat_fun.h	67.64%	11 Missing ⚠️
...tTests/material_model_tests/test_material_common.h	56.25%	7 Missing ⚠️
Code/Source/solver/post.cpp	81.25%	6 Missing ⚠️
Code/Source/solver/mat_models.h	0.00%	4 Missing ⚠️
Code/Source/solver/fsi.cpp	50.00%	3 Missing ⚠️
Code/Source/solver/mat_models.cpp	92.10%	3 Missing ⚠️
Code/Source/solver/nn_elem_props.h	0.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #640      +/-   ##
==========================================
- Coverage   73.41%   72.88%   -0.53%     
==========================================
  Files         269      272       +3     
  Lines       40337    40148     -189     
  Branches     6755     6718      -37     
==========================================
- Hits        29612    29263     -349     
- Misses      10482    10663     +181     
+ Partials      243      222      -21     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[aabrown100-git]

@dseyler the viscous stress performance improvements look good, but why is the wedge element shape function included in this PR? Also, is it possible to add a wedge element test case?

[dseyler]

@dseyler the viscous stress performance improvements look good, but why is the wedge element shape function included in this PR? Also, is it possible to add a wedge element test case?

@aabrown100-git Maybe I should have opened a 2nd PR for that, but it's just a one-line change that @ktbolt noted was mistranslated from Fortran. The connection between the two PRs is that the viscosity model now only computes once for the first gauss point of each linear element as the viscosity arrays are the same for all Gauss points in the element. This would have given incorrect results for wedge elements which had the wrong lShpF flag. Several other places in the code also skip unnecessary per-gauss point computation for linear elements and would have given incorrect results for wedge elements as well

Definitely a good idea to add some wedge element test cases in the future, as there are currently none. Do you know of anyone working with them?

[aabrown100-git]

I thought @msbazzi was working with them?

If you can come up with a test case, I would create a new PR with your bug fix and a test case.

[ktbolt]

A user on the SV Forum just sent me a mesh that uses wedges for the boundary layer. svMultiPhysics fails reading the mesh though; I will investigate.

[dseyler]

@aabrown100-git @ktbolt Are you suggesting creating a wedge element test case to add to this PR, or reverting the bug fix and pushing this PR through without it, which would cause the viscosity model to produce incorrect results for wedge elements (as 30 other places in the codebase already do)?

I haven't worked with wedge elements before, so would defer to someone else if they already have a good wedge element case to add

[ktbolt]

@dseyler Go ahead and merge this PR. I will create a separate Issue for adding a wedge element CI test and for fixing the problem reading in wedges if it is indeed a bug and not something wrong with the user's file.

[aabrown100-git]

Approving!

[ktbolt]

@aabrown100-git @dseyler Note that you can created wedges using vtkLinearExtrusionFilter.

[lassoan]

@aabrown100-git The new "CFD mesh generator" module in VMTK Slicer extension can generate wedges as boundary layers (with multiple iterations of extrusion, smoothing, detangling). This can be useful if you want to test on complex, realistic meshes.

[dseyler]

Addressing comments here and the many changes that have been made to the branch:

• struct3d/2d and the ustruct per-Gauss-point assembly routines have been refactored using Eigen. This considerably condensed much of the arithmetic and made it much more readable. To facilitate this, the following changes were made:
  • eigen_view helper functions were added to construct read-only eigen views, reducing the amount of Eigen jargon in the assembly routine.
  • Aliases for fixed-size and bounded dynamic Eigen matrices were added in mat_fun
  • compute_pk2cc is now templated on nsd to accept fixed-size Eigen. An array overload still exists, used only in post.cpp. If preferred, this overload could be removed and the arrays can be converted to Eigen in post.cpp with a few loc before being passed into the templated functions.
  • const int MaxNoN = 27 is added in consts.h and used to bound dynamic Eigen matrices. @zasexton confirmed this is the appropriate location. This value was previously defined in several different files and is now consolidated in one place. If larger element sizes are supported in the future, this will need to be updated, but it is now much easier to catch than before. In the future, this could be linked to nn_elem_props.h so that it automatically updates if larger elements are added.
  • recompute_visc is checked in the assembly routine and not passed to the viscosity function.
  • Scratch arrays for Kvis_u and Kvis_v were previously declared too far upstream of their use. Rather than move this into the struct_3d/2d where complicated resizing logic would be necessary to make use of recompute_visc, they are now owned by a new ViscousReponse class which hides the resize from the assembly routine. The viscosity call sites now read visc.update(..) and Svis/Kvis_u/Kvis_v are accessed as visc.S()/visc.du/visc.dv. Maybe this is a bit more of a refactor than necessary, but I think is much more readable.
  • Faster tensor operations were added in mat_fun.h. Maybe this should be a separate PR, but they were already in my branch that I merged with the viscosity branch to address @ktbolt and @michelebucelli's request to move Eigen as far up as possible.

A few more notes:

• Chasing Eigen conversions further upstream gets complicated, as it would require editing nn::gnn which has a wide impact on many modules. I also think it would have diminishing returns, as the functions edited in this branch are called per-Gauss point while the rest of the assembly routine is per-element or less frequent.
• Many of these same changes could be made to b_struct* and would make those functions easier to read, but not necessarily much faster.
• While implementing Eigen in ustruct, I found many bugs in ustruct_2d where arrays are clearly indexed out of bounds, either due to incorrect indexing conversion from FORTRAN or from copying 3d code into the 2d routines. There are currently no 2d ustruct test cases, so these bugs would not have been caught. I tried to fix all of the obvious bugs during the Eigen conversion but can not promise that more don't exist.

[zasexton]

General question, when building the solver and running the Eigen Matrices and Tensors are you compiling with the -march flag in cmake? I believe this may be necessary for the Eigen vectorization optimizations.

[zasexton]

It would be nice if the maxNoN was of type size_t since we would otherwise expect this to always be a positive integer value; however I understand that Eigen templating is built around int type for dimension sizing and because dynamic runtime allocation is represented by -1.

[dseyler]

General question, when building the solver and running the Eigen Matrices and Tensors are you compiling with the -march flag in cmake? I believe this may be necessary for the Eigen vectorization optimizations.

I was not. But it looks like Apple Silicon vectorizes with Eigen by default with ARM NEON. Is there a reason why this isn't the default on Sherlock though?

[dseyler]

It would be nice if the maxNoN was of type size_t since we would otherwise expect this to always be a positive integer value; however I understand that Eigen templating is built around int type for dimension sizing and because dynamic runtime allocation is represented by -1.

Sounds good. A lot of the templates were previously on size_t and I switched them to int per Claude's recommendation (for the reason you mentioned). I can switch them back though. I think it should be safe to use size_t as long as it's always positive. Will this be fine for the places where the templated parameter has arithmetic in it, i.e. Matrix<3*(nsd-1)>?

[zasexton]

General question, when building the solver and running the Eigen Matrices and Tensors are you compiling with the -march flag in cmake? I believe this may be necessary for the Eigen vectorization optimizations.

I was not. But it looks like Apple Silicon vectorizes with Eigen by default with ARM NEON. Is there a reason why this isn't the default on Sherlock though?

It seems like someone made a note of this at some point

svMultiPhysics/CMakeLists.txt

Line 48 in c9838bd

# Add -march=native to CMAKE_C_FLAGS and CMAKE_CXX_FLAGS for Eigen performance

In general, I'd suspect that the -march flag is not used by default because it specifies hardware specific accelerators on the machine it is compiled on. This might cause problems if you compile against a specific CPU architecture but then you run a simulation across nodes that might not share the exact same type of architecture. I think we currently have two different CPU architectures on the sherlock partition (MLN and SKX).

[zasexton]

It would be nice if the maxNoN was of type size_t since we would otherwise expect this to always be a positive integer value; however I understand that Eigen templating is built around int type for dimension sizing and because dynamic runtime allocation is represented by -1.

Sounds good. A lot of the templates were previously on size_t and I switched them to int per Claude's recommendation (for the reason you mentioned). I can switch them back though. I think it should be safe to use size_t as long as it's always positive. Will this be fine for the places where the templated parameter has arithmetic in it, i.e. Matrix<3*(nsd-1)>?

So long as any arithmetic is properly guarded for underflow then yes I believe that size_t should work. But before making any changes we can see what @ktbolt and @michelebucelli think about the nsd type preferences for maxNoN?

[michelebucelli]

So long as any arithmetic is properly guarded for underflow then yes I believe that size_t should work. But before making any changes we can see what @ktbolt and @michelebucelli think about the nsd type preferences for maxNoN?

@zasexton @dseyler I generally lean in favor of using unsigned types (like size_t) for quantities that are supposed to be non-negative. I see this first and foremost as a declarative practice, as in: the person writing the code declares to the person reading it that a number is assumed to be non-negative, with the variable type (as opposed to a documentation comment) documenting this intent.

I think this applies to maxNoN and nsd too, so I would prefer size_t over int for both of them.

Underflow is a possible concern, but given the semantics of these variables it should be fine without guards most of the time, and can be guarded with static_asserts in those cases where it isn't (e.g. before declaring something of size nsd - 2, one could place a static_assert(nsd >= 2, ...) just to be safe).

The guard is used to occasionally enforce a stricter constraint than just non-negativity (e.g. it must be greater than 2 in previous example), and this is better than using ints, which amounts to no compile-time constraint at all (or would require static_assert constraints in far more places).

[michelebucelli]

Thanks @dseyler! I did another review pass.

I think using Eigen the code looks much cleaner and more readable!

I still have some doubts about the current implementation of viscous stress caching across Gauss points, sorry if I keep returning to that (and for occasionally contradicting my own previous suggestions).

Chasing Eigen conversions further upstream gets complicated, as it would require editing nn::gnn which has a wide impact on many modules. I also think it would have diminishing returns, as the functions edited in this branch are called per-Gauss point while the rest of the assembly routine is per-element or less frequent.

I think long term it might be a good idea to use Eigen for every small, dense linear algebra task (parallel and/or sparse linear algebra is a different question entirely). I can see your point about diminishing returns in terms of efficiency, but I think the returns would not be diminishing in terms of code hygiene and maintenance burden. I agree that such a large refactor would be out of scope for this PR, so I'm fine with what you currently did (and also I appreciate how this PR sets a good example of the improvements that can be obtained).

I think we might make a bit of a collective effort to slowly migrate to Eigen the parts of code that we touch when we get the opportunity (and maybe remember to gently suggest to do this in otherwise unrelated PR reviews 😅).

[dseyler]

I still have some doubts about the current implementation of viscous stress caching across Gauss points, sorry if I keep returning to that (and for occasionally contradicting my own previous suggestions).

Thanks for reviewing all of this for a second time! I agree that the viscous stress caching is a bit cumbersome and am not married to it. The ultimate goal of that feature was to save people time by making our simulations faster, but if it costs us more time to review and rewrite than it saves, then it might defeat the purpose 😅

[michelebucelli]

Approved, apart from one final minor renaming suggestion.