sandialabs/LCM

Rework LCM to not use deprecated STK code

ikalash opened this issue ยท 12 comments

See #85 for more detail. LCM needs to be reworked to not use deprecated code, in a similar way as was done for Albany master here: sandialabs/Albany@9d593cd in early June.

@ikalash You're right that the changes should be very similar to the Albany commit you referenced. Absolutely everything that is going away has a deprecation warning (mostly compile-time, but things like the Field template parameter tags are now run-time), so that can guide you on what to remove. The files CoordinateSystems.hpp, TopologyDimensions.hpp, and FieldTraits.hpp will eventually be deleted, so their includes should be removed. Also, you should include a call to either MetaData::use_simple_fields() or StkMeshIoBroker::use_simple_fields() soon after their construction to make sure that any Fields that are auto-registered on your behalf (like coordinates) during mesh read conform to the new style. This will also change any usage of the old-style Fields into a run-time error and prevent accidental regressions before the legacy code support is removed from STK. If you take care of all of this, then you shouldn't notice at all when we remove the legacy Field support later this year.

Feel free to ping me here if you have any questions at all. I'll be happy to help, especially since this conversion is somewhat of a large burden to place on our users.

Thank you @djglaze, this is very helpful! I will likely take you up on your offer to answer questions about the deprecations. I am getting ready to attend a foreign conference next week and will be out for the next 2 weeks, but I plan to look at this after I get back.

@lxmota : this deprecated code has been removed, and this is why the LCM nightlies are failing. I will try to fix this soon.

@ikalash Ah ok. Since you have done this before that would be great, but let me know if I can help; or if you tell me what I need to fix I'll do it.

@lxmota : I actually haven't done it before, but Mauro has and pointed me to the PR where the changes were done.

@ikalash Ok, let me know how I can help.

I have finally started working on this, and had a question for @bartgol and @djglaze . While the changes to STK are pretty clear (getting rid of VectorField, TensorField layouts, removing stk::mesh::Cartesian, etc.), I am a bit worried that some aspects of LCM's design are not amenable to the STK changes. It looks like back in 2021, @bartgol you refactored the STK discretization classes in Albany to get rid of BucketArrays, replacing it with code in Albany_STKUtils.hpp. This was not done in LCM, so we still have the BucketArray stuff, overloaded from Shards: https://github.com/sandialabs/LCM/blob/main/src/disc/stk/Albany_BucketArray.hpp . My question is, can BucketArray be made to work with the STK changes? There is code that makes use of stk::mesh::Cartesian, for example, which would not work anymore (https://github.com/sandialabs/LCM/blob/main/src/disc/stk/Albany_BucketArray.hpp#L90). I've started trying to convert LCM to use Albany_STKUtils.hpp but haven't finished so I might hit some snags. Is there any reason why it would not be possible to modify LCM to use Albany_STKUtils.hpp and get rid of buckets? One difference in LCM is that we have Cauchy stress fields which are the Tensor layout, so it needs to work correctly. I believe the refactor to take out buckets in Albany was independent of changes to refactor Albany to have block discretizations.

Looking at it more, to switch away from using buckets, it seems some fairly substantial changes are required, such as passing in dof_manager in place of buckets to fillVector in Albany::STKContainerHelper (see https://github.com/sandialabs/Albany/blob/master/src/disc/stk/Albany_STKFieldContainerHelper.cpp#L28 in Albany vs. https://github.com/sandialabs/LCM/blob/main/src/disc/stk/Albany_STKFieldContainerHelper_Def.hpp#L82 in LCM), and possibly other files that call it. If we can avoid this by keeping the buckets, I would maybe prefer that, but it would be good to get some insight into whether things can be made to work with buckets or not.

If you still want to stuff data in BucketArray (which is just a shards array), I think you can still do it, no? I suppose you can use any shards tag in place of stk::mesh::Cartesian. E.g., any of the tags in PHAL_Dimension.hpp would work, no? Maybe I'm not understanding the core issue here...

@bartgol : thanks for replying and sorry for not being clear. I haven't worked with shards/BucketArray in a very long time, and when I noticed all the specializations within it for ScalarFieldType, VectorFieldType, etc. (see https://github.com/sandialabs/LCM/blob/main/src/disc/stk/Albany_STKFieldContainerHelper_Def.hpp#L16-L68) and also usage of stk::mesh::Cartesian (https://github.com/sandialabs/LCM/blob/main/src/disc/stk/Albany_STKDiscretization.cpp#L1728-L1730) when assigning data (https://github.com/sandialabs/LCM/blob/main/src/disc/stk/Albany_STKDiscretization.cpp#L1728-L1730), I got concerned about the approach being viable.

I guess in Albany_STKFieldContainerHelper_Def.hpp I just need to have a way to tell the code whether a 2D or 1D array is being provided to access it properly, right? This wouldn't come from the type of argument anymore, which will always be based on STKFieldType rather than ScalarFieldType, VectorFieldType, etc. now.

In Albany_STKDiscretization.cpp, I guess I need to modify array.assign<ElemTag, NodeTag, CompTag> and other calls involving ElemTag,NodeTag, etc. that are not of type stk::mesh::Cartesian. What would be the appropriate tags to use here?

I guess in Albany_STKFieldContainerHelper_Def.hpp I just need to have a way to tell the code whether a 2D or 1D array is being provided to access it properly, right? This wouldn't come from the type of argument anymore, which will always be based on STKFieldType rather than ScalarFieldType, VectorFieldType, etc. now.

Yes, that's correct. The type of the input field is no longer different, so you need a different way to query the rank. You can use stk::mesh::field_scalars_per_entity (in one of its flavors). Given that, and knowing how many QPs are in an elem, you can infer if it's a scalar or vector field (that is, assuming you never have vector fields with vector-extent 1...). It gets tricky to distinguish between a 9dim vector and a 3x3 tensor: for stk, these are the same. But I think that should not be an issue, since you'd store both vec-9 and tens-3x3 as the same kind of field in the stk mesh, so no real issue when accessing the stk data.

In Albany_STKDiscretization.cpp, I guess I need to modify array.assign<ElemTag, NodeTag, CompTag> and other calls involving ElemTag,NodeTag, etc. that are not of type stk::mesh::Cartesian. What would be the appropriate tags to use here?

I would use the tags we define in PHAL_Dimension.hpp. They seem like the natural choice.

Thanks @bartgol , this is very helpful! I think what I will do is try to rework LCM to use stk::mesh::field_scalars_per_entity and PHAL_Dimension.hpp objects, and make sure everything works with that first, before switching over to STKFieldType. This should minimize the number of pieces that are moving at the same time, and error, hopefully. I might ping you by email if I have further questions / hit snags, since you are likely more up-to-date with the relevant stk routines that can be helpful than I am.