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.08% 72.88% -0.21%
==========================================
Files 268 273 +5
Lines 40270 40172 -98
Branches 6738 6722 -16
==========================================
- Hits 29433 29281 -152
- Misses 10594 10670 +76
+ Partials 243 221 -22 ☔ 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. |
…me Eigen expressions
|
Addressing comments here and the many changes that have been made to the branch:
A few more notes:
|
|
General question, when building the solver and running the Eigen Matrices and Tensors are you compiling with the |
|
It would be nice if the |
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? |
Sounds good. A lot of the templates were previously on |
It seems like someone made a note of this at some point Line 48 in c9838bd In general, I'd suspect that the |
So long as any arithmetic is properly guarded for underflow then yes I believe that |
@zasexton @dseyler I generally lean in favor of using I think this applies to 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 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 |
michelebucelli
left a comment
michelebucelli
left a comment
There was a problem hiding this comment.
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 😅).
|
|
||
| const int maxNSD = 3; | ||
|
|
||
| const int maxNoN = 27; // Max node count in nn_elem_props.h |
There was a problem hiding this comment.
-
I imagine you have discussed this internally already, but: why not place this in
nn_elem_props.h, since that's where the information that this number is derived from is also stored? -
Also, how does this relate to the FE basis refactoring recently implemented by @zasexton? Do you envision there will be a place in the new infrastructure to save/query this information?
-
It might be a good idea to declare this (and maybe other variables in this file) as
constexprinstead ofconst.It shouldn't make a difference most of the time, since I expect the compiler will be able to figure it out by itself, but it makes it more explicit to the reader that these are supposed to be used in compile-time expressions (e.g. template as arguments).
| // The element routines know their dimension at compile time and call the | ||
| // template directly, so instantiate the dimensions the solver supports. |
There was a problem hiding this comment.
Minor: I suggest dropping the comment, because it reads a bit abstract and AI-ish without carrying much meaning. If we want to keep it, I would suggest rephrasing to something more concrete, e.g:
// Explicitly instantiate compute_pk2cc for 2D and 3D.| fl_3D(1, i) = fl(1, i); | ||
| fl_3D(2, i) = fl(2, i); | ||
| } | ||
| auto F_3D = mat_fun::convert_to_eigen_matrix<Matrix<3>>(F); |
There was a problem hiding this comment.
I suggest making this const (the compiler probably figures it out on its own, but it might be better to be explicit).
| /// @brief Writable Eigen view of a whole Array, sharing its storage. | ||
| inline Eigen::Map<Eigen::MatrixXd> | ||
| eigen_view_mut(Array<double>& A) { | ||
| return {A.data(), A.nrows(), A.ncols()}; | ||
| } |
There was a problem hiding this comment.
Minor: I suggest the more explicit name eigen_view_mutable.
| template<size_t nsd> | ||
| template<int nsd> | ||
| using Matrix = Eigen::Matrix<double, nsd, nsd>; |
There was a problem hiding this comment.
I'm not sure myself about this suggestion, but I think there's an argument to be made for this to be called Tensor instead of matrix (and the current Tensor alias below be renamed to Tensor4).
This type is used to represent square matrices with as many elements as there are spatial dimensions (so the name should maybe be more restrictive than just Matrix to reflect this). In the constitutive modelling code this is used to represent objects that, mathematically, are tensors rather than general square matrices. I think changing the name would make this more apparent.
| template<size_t nsd> | ||
| template<int nsd> | ||
| using Matrix = Eigen::Matrix<double, nsd, nsd>; | ||
|
|
||
| template<size_t nsd> | ||
| template<int nsd> | ||
| using Tensor = Eigen::TensorFixedSize<double, Eigen::Sizes<nsd, nsd, nsd, nsd>>; | ||
|
|
There was a problem hiding this comment.
I suggest to give Doxygen documentation to these aliases (especially with reference to their intended use and whether they are fixed- or dynamic-sized).
| /** | ||
| * @brief Contracts two 4th order tensors A and B over two dimensions, | ||
| * | ||
| * @brief Contracts two 4th order tensors A and B over two dimensions. |
There was a problem hiding this comment.
I think it would be helpful to report the expression on the contraction (e.g. in Einstein notation) in the documentation, to make sure this is clear.
| /// @brief Array-based overload, for callers whose dimension is a run-time value. | ||
| void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, | ||
| const dmnType &lDmn, const Array<double> &F, const int nfd, | ||
| const Array<double> &fl, const double ya_f, | ||
| const double ya_s, const double ya_n, Array<double> &S, | ||
| Array<double> &Dm, double &Ja); |
There was a problem hiding this comment.
Are there any call sites left where this is needed? If so, would it be possible to have them use Matrix too, and drop this overload altogether? Or is this needed to dynamically dispatch at runtime to the correct dimension (since com_mod.nsd is a run-time variable)?
The less overloads there are, the easier it is for the reader to figure out which function is being called (especially for functions with this many arguments).
There was a problem hiding this comment.
It's used in post.cpp and the unit tests. The array overload could be removed and just dispatched to the templated version in post.cpp (with a few lines converting arrays to Eigen before calling compute_pk2cc)
| /// @tparam nsd Number of spatial dimensions. | ||
| /// | ||
| template <int nsd> | ||
| class ViscousResponse { | ||
| public: | ||
| /// @brief Evaluate the domain's viscosity model at this Gauss point. | ||
| /// | ||
| /// @param[in] lDmn Domain, supplying the viscosity model and its parameters. | ||
| /// @param[in] eNoN Number of element nodes. | ||
| /// @param[in] Nx Shape function spatial derivatives. | ||
| /// @param[in] vx Velocity gradient. | ||
| /// @param[in] F Deformation gradient. | ||
| void update(const dmnType& lDmn, const int eNoN, const Array<double>& Nx, | ||
| const Matrix<nsd>& vx, const Matrix<nsd>& F); | ||
|
|
||
| /// @brief Viscous 2nd Piola-Kirchhoff stress. | ||
| const Matrix<nsd>& S() const { return Svis_; } | ||
|
|
||
| /// @brief Tangent w.r.t. displacement. du(i*nsd + j, a, b) is the (i,j) | ||
| /// entry of the block coupling nodes a and b. | ||
| double du(const int ij, const int a, const int b) const { return Kvis_u_(ij, a, b); } | ||
|
|
||
| /// @brief Tangent w.r.t. velocity. | ||
| double dv(const int ij, const int a, const int b) const { return Kvis_v_(ij, a, b); } | ||
|
|
||
| private: | ||
| Matrix<nsd> Svis_; | ||
| Array3<double> Kvis_u_, Kvis_v_; | ||
| }; |
There was a problem hiding this comment.
If I understand this right, this class is meant to be a helper that assembly routines use to walk over the mesh elements and evaluate a certain quantity (viscous stress) at their integration points. Is this correct?
If it is, I think this has some analogy with the ActiveStress::Evaluator class I have implemented in #650, and I expect a similar pattern to pop up in a number of other places.
I think it might be a good idea to abstract this away and establish a common infrastructure for "element walkers"/"integration point evaluators", to avoid repeated work, ensure consistency, and facilitating adding new instances of this pattern.
I propose we discuss this after this PR and #650 are merged.
There was a problem hiding this comment.
I mostly implemented this to hide the whole Kvis_u/v resizing logic from struct3d/2d, but I agree it would be nice to have some common infrastructure for this kind of thing. Ideally, I feel like the different contributions to the tangent and residual should be able to be abstracted away so that the whole assembly routine can be separated into viscous, constitutive, prestress, etc. contributions instead of a long routine of disordered arithmetic.
Right now construct_d_solid loops over each element and each Gauss point and then runs struct3d to assemble lK and lR per gauss point. I feel like a lot of this code would be simpler if construct_d_solid (or a future general assembly function) just looped over the elements, and then physics-specific assembly kernels looped over Gauss points. This would make it much easier to separate quantities that need to be computed once per element (viscosity, shape function derivatives for linear elements) and those that need to be computed per Gauss point and would remove the need for the recompute_visc flag
| // Viscous 2nd Piola-Kirchhoff stress and tangent contributions. | ||
| // Reuse from the previous Gauss point when shape function gradients | ||
| // are constant within an element (e.g. linear triangles, tetrahedra). | ||
| static mat_models::ViscousResponse<2> visc; | ||
| if (recompute_visc) { | ||
| visc.update(dmn, eNoN, Nx, vx, F); | ||
| } |
There was a problem hiding this comment.
I think this works under the assumption that repeated calls to this method happen in a particular order, which is currently true, but that this method itself has no control over. Should this change in the future, it might give rise to bugs that are very hard to find (and even to notice, probably).
A possible better design could be one where the methods struct_2d and struct_3d are not responsible for evaluating/updating the viscous stress, but receive a ViscousResponse object in input, and this object is managed by the caller. This way, the same function that loops over the elements is the one that decides whether to recompute viscosity or not.
Another possible design would be one where this is delegated to ViscousResponse itself, i.e. its update method decides whether it really has to update or it can silently reuse the already computed quantities. This however would require passing more information into the ViscousResponse class, and I'm not sure we want to do that. It's something we can discuss as part of the proposed evaluator abstraction that I'm mentioning in some other comment.
Similar considerations apply to other places where this class is used (also in ustruct).
(Sorry for going back and forth on this issue, I am struggling a bit to find a clean and safe solution that still allows reuse of the viscous tensors).
There was a problem hiding this comment.
Didn't read this before my response to your comment above on the ViscousResponse class.
But if I understand right, we both agree a lot of this complexity comes from the fact that the boundary between constructd_solid and struct_2/3d is at the level of the Gauss points rather than the element.
If the boundary were in the element loop, the update method could decide what to do depending on the element type, but currently it doesn't have that info since the element loop lies one layer up in constructd_solid. If visc.update were called once per element, the method would not need to have a recompute_visc flag passed in
There was a problem hiding this comment.
This is a good way of framing the issue. However, even if visc.update were called once per element, it would need to decide whether to compute a different viscous tensor for every integration point or not, which would depend on information that the ViscousResponse currently doesn't have.
Because of this, I still think the two problems (whether ViscousResponse works per element or per Gauss node, and how to avoid repeated computations) are not exactly the same.
I think for the sake of this PR the first solution I propose (evaluating viscosity in the caller of struct_2d/3d) would be a reasonable solution, but we may want to revisit this if we abstract the integration point evaluator pattern away (which I think we should, but after this).
There was a problem hiding this comment.
Sounds good! I'll go with your first solution
There was a problem hiding this comment.
Since ViscousResponse is templated on runtime nsd, if it's moved into the caller, then either 1 (easiest): both need to be declared in construct_d_solid, construct_u_solid and construct_d_fsi or 2: those functions could also be templated on nsd as well.
Easiest:
/// Declare both. Use only one depending on whether struct3d or struct 2d is called
mat_models::ViscousResponse<2> visc2;
mat_models::ViscousResponse<3> visc3;
Ideally, this would be cleaner. Do you have thoughts?
There was a problem hiding this comment.
Looked into this a bit more:
Computing the ViscousResponse before struct_3d/2d and passing it in as an argument would require hoisting the computation of its arguments F and vx out of struct3d/2d as well.
I think this is a good plan to do eventually, since F and vx are also uniform over linear elements (they only depend on shape function derivatives), so it's a little odd that viscosity is computed once per element but F and vx are computed for each Gauss point even though F and vx are arguments of visc.update.
If this were implemented, then the recompute_visc check would be entirely contained in construct_dsolid() and read:
for e < lM.nEL {
for g < lM.nG {
if nsd == 2:
if (g == 0 || !lM.lShpF) {
nn::gnn(Nx)
Nxm = eigen_view<2>(Nx)
vx = ...
F = ...
visc.update()
struct2d(...,visc)
if nsd == 3
.
.
.
I think at that point it would make sense to convert gnn to Eigen and template it on nsd so the vx and F computation is a little more legible. gnn was also one of the biggest consumers of runtime from my earlier profiling work, but I didn't look into it because it touches so many other modules. Also, just looking at the function, it is a LOT of unrolled matrix multiplications.
How does this sound for a plan:
-
Either pass recompute_visc as a parameter into visc.update() (Your suggestion 2) or scrap the recompute_visc idea for now and come back to it in a later PR. I don't feel strongly either way.
-
At a later date, convert
gnnto Eigen and movevisc.updateout ofstruct3d.We could create similar methods forcompute_pk2cc's passive and active contributions to S, i.e. passive.update() and active.update() as well, using the general Evaluator class you suggested.
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 😅 |
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