Templated viscous stress functions + linear shape function fix - #640
Templated viscous stress functions + linear shape function fix #640dseyler wants to merge 11 commits into
Conversation
…ted to false for wedge elements
![claude[bot]](https://avatars.githubusercontent.com/in/1236702?s=60&v=4){kind=small}
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #640 +/- ##
==========================================
- Coverage 73.14% 73.08% -0.07%
==========================================
Files 259 259
Lines 39516 39514 -2
Branches 6673 6667 -6
==========================================
- Hits 28904 28877 -27
- Misses 10368 10393 +25
Partials 244 244 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
@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 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? |
|
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. |
|
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. |
|
@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 |
|
@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 @dseyler Note that you can created wedges using |
|
@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. |
There was a problem hiding this comment.
🟡 Changes recommended
The updated viscosity kernels introduce silent-failure paths and reduced singular-matrix error handling that can leave outputs stale or produce NaNs in edge cases.
Once you've addressed the issues Copilot identified, you can request another Copilot review.
Pull request overview
This PR improves structural-mechanics assembly performance by templating viscous stress kernels on nsd and enabling reuse of viscous stress/tangent terms across Gauss points when shape-function gradients are constant, while also fixing an element-property bug for wedge elements (WDG) where gradients were incorrectly treated as constant.
Changes:
- Refactors Newtonian/potential viscosity models in
mat_models.cpptonsd-templated implementations using Eigen stack-based views. - Plumbs
Svis,Kvis_u,Kvis_v, and arecompute_viscflag throughstruct_*,ustruct_*, and FSI assembly loops to reuse viscosity computations for linear elements. - Fixes wedge element properties by setting
lShpF=falseso shape gradients/Jacobians are recomputed per Gauss point.
File summaries
| File | Description |
|---|---|
| Code/Source/solver/ustruct.h | Extends ustruct_*_m APIs to accept viscous stress/tangent buffers plus recompute_visc. |
| Code/Source/solver/ustruct.cpp | Allocates/reuses viscosity buffers across Gauss points and passes recompute_visc into viscosity computation. |
| Code/Source/solver/sv_struct.h | Extends struct_2d/struct_3d APIs to accept viscosity buffers plus recompute_visc. |
| Code/Source/solver/sv_struct.cpp | Allocates/reuses viscosity buffers across Gauss points and passes recompute_visc into viscosity computation. |
| Code/Source/solver/nn_elem_props.h | Sets WDG lShpF=false to avoid incorrectly reusing shape gradients within wedge elements. |
| Code/Source/solver/mat_models.h | Updates compute_visc_stress_and_tangent signature to include recompute_visc and improves its API doc. |
| Code/Source/solver/mat_models.cpp | Implements nsd-templated viscosity kernels and conditional recomputation/dispatch by nsd. |
| Code/Source/solver/mat_fun.h | Adds fixed-size Eigen overloads for mat_symm and mat_dev. |
| Code/Source/solver/fsi.cpp | Allocates/reuses viscosity buffers and forwards them into solid element routines with recompute_visc. |
Review details
Suppressed comments (1)
Code/Source/solver/mat_models.cpp:1759
- Same dispatch issue as the Newtonian branch: if
F.nrows()is not 2 or 3, the Potential case becomes a no-op and can leaveSvis/Kvis_*stale.
if (F.nrows() == 3) {
compute_visc_stress_potential<3>(lDmn.solid_visc.mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v);
} else if (F.nrows() == 2) {
compute_visc_stress_potential<2>(lDmn.solid_visc.mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v);
}
- Files reviewed: 9/9 changed files
- Comments generated: 3
- Review effort level: Lite
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
michelebucelli
left a comment
michelebucelli
left a comment
There was a problem hiding this comment.
Thank you @dseyler! I've left some suggestions.
As a general comment, I appreciate that Eigen brings some nice performance improvements, and that some of these changes also add some efficiency by avoiding unnecessary repeated allocations. However, there's a risk that both these kinds of changes increase the amount of "noise" in the code, by which I mean that the functions not only do the mathematical computation that they are supposed to do, but also spend quite a few lines in data format conversion and general bookkeeping, which may obfuscate their meaning a bit.
Some of the comments below go into the direction of mitigating this kind of issue (they are of course very much open to debate). I think if we can figure out some nice solution to this, then it can be used in the other performance improvements that you've worked on.
| } | ||
|
|
||
| if (g == 0 || !fs_1[0].lShpF) { | ||
| // Viscosity is constant at all Gauss points for linear elements |
There was a problem hiding this comment.
Minor suggestion:
| // Viscosity is constant at all Gauss points for linear elements | |
| // Viscosity is constant at all Gauss points for linear tetrahedral elements |
Linear hexahedral elements are actually bi- or trilinear, so I don't think that this would apply to them, and I think it's a good idea to leave a reminder in the comment.
(I know that technically those are not linear elements, but they are sometime called that, so it might be better to be redundantly explicit just in case).
| * Fixed-size overload for the Eigen matrices used by the element kernels. | ||
| * | ||
| * @tparam nsd Number of spatial dimensions. | ||
| * @param[in] A Second order tensor. |
There was a problem hiding this comment.
Very minor:
| * @param[in] A Second order tensor. | |
| * @param[in] A second order tensor. |
Same for the function below.
| /// @brief Largest element node count the fixed-size views below allow (HEX27). | ||
| constexpr int MAX_ELEMENT_NODES = 27; |
There was a problem hiding this comment.
I'm not sure I have a clean solution to this, but: this number depends on the "largest" local finite element basis supported by the library, and this information is not logically pertinent to this module (the material model evaluation) but to the finite element basis module.
I think having this constant here introduces possible unexpected bugs (e.g. if larger finite element spaces are added in the future), and in general introduces implicit coupling between the two modules (implicit in the sense that the coupling is somewhat hidden).
One somewhat cleaner solution to this, I think, would be to move this definition to the basis function module (source/Code/FE/Basis). I am not sure what the appropriate file for that would be, but maybe @zasexton has a clearer picture.
| Eigen::Map<const Matrix<nsd>> F_map(F.data()); | ||
| Eigen::Map<const Matrix<nsd>> vx_map(vx.data()); | ||
| Eigen::Map<const NodalMatrix<nsd>> Nx_map(Nx.data(), nsd, eNoN); |
There was a problem hiding this comment.
If I understand this right, these lines are constructing Eigen views over the data stored by the arrays, so that the subsequent operations use Eigen for efficiency but don't do any copy.
If that is right,
- Would it make sense to change the code in the caller to directly use
Eigenstructures for the local data, instead of customArrays? This would be a somewhat larger change perhaps, but I think it might make sense, reduce the (presumably very small, if even detectable) overhead associated to the construction of these views, and make the code a bit simpler (this function would just compute the viscosity, without any boilerplate code related to the technicalities of conversion). - Should you choose not to follow the above suggestion, I would recommend writing a small helper function to do this, as a form of synctactic sugar, to turn these lines into e.g.
auto F_eigen = to_eigen_map(F).
A similar suggestion applies to compute_visc_stress_newton.
| template <int nsd> | ||
| void compute_visc_stress_potential(const double mu, const int eNoN, const Array<double>& Nx, | ||
| const Array<double>& vx, const Array<double>& F, | ||
| Array<double>& Svis, Array3<double>& Kvis_u, Array3<double>& Kvis_v) { |
There was a problem hiding this comment.
If I interpret this right, you removed the definition of this function and of compute_visc_stress_newton from the header because they are not meant to be called from outside this file.
If that is correct, I think it is a good decision, but I would also place both functions in the anonymous namespace above (my understanding of anonymous namespaces is that, among other things, they're a way of specifying that a certain name should only have internal visibility).
| if (!recompute_visc) { | ||
| return; | ||
| } |
There was a problem hiding this comment.
I think that it would be better if this check was performed by the caller of this function, rather than by the function itself.
In other words, I think this function should always do what its name says (compute the viscosity), and the caller, who is aware of the context, should be the one to decide whether the viscosity needs to be recomputed (and thus the function needs to be called) or not.
| Array<double> Svis(nsd,nsd); | ||
| Array3<double> Kvis_u(nsd*nsd,eNoN,eNoN), Kvis_v(nsd*nsd,eNoN,eNoN); |
There was a problem hiding this comment.
It seems to me that these arrays are not used by this function itself, but rather they are only used as "scratch" writing space for struct_3d and struct_2d, and they are declared here to avoid repeated allocations.
One downside is that this forces us to declare variables in a scope where those variables are not meaningful, which might make it harder to interpret the code.
An alternative solution to that would be to declare these variables inside the functions that need them (struct_3d and struct_2d) but declare them as static. This way, they will not be reallocated between subsequent calls, but they will remain local to the function's scope. You might have to write a bit of logic to ensure that they always have the correct size between calls, but I think it could be a good solution.
Current situation
Address issues #633 and #634, which involve templating viscosity models on
nsdand 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_potentialare now templated onnsdinmat_models.cppcompute_visc_stress_and_tangentdispatches to templated models bynsdmat_symmandmat_devare alsonsdtemplated inmat_fun.hSvis,Kvis_u, andKvis_vare 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_viscflag is added to indicate when viscosity computation can be skipped (shape function gradient are uniform for all Gauss points within an element, soKvis_u,Kvis_v, andSvisvalues can be reused).compute_visc_stress_and_tangentskips viscosity computation whenrecompute_visc == Falsecompute_visc_stress_and_tangentso viscosity arrays are zeroed and do not carry the previous element's values when an element doesn't have a viscosity model definedstruct_2d,struct_3d,ustruct_2d_m,ustruct_3d_mnow require four new arguments:Svis,Kvis_u,Kvis_v, andrecompute_visccompute_visc_stress_and_tangentnow takesrecompute_visclShpFshould 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