Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[MPM] Moving enforcement of nonconforming SLIP from MPMBoundaryRotationUtility to to the relevant Condition #12018

Merged
merged 24 commits into from
Aug 9, 2024

Conversation

gzjing
Copy link
Contributor

@gzjing gzjing commented Feb 5, 2024

Currently, the MPMBoundaryRotationUtility is used to enforce both conforming and non-conforming slip Dirichlet conditions. While this makes sense in the conforming case, for which the contributions of all conditions and elements need to be modified, this makes little sense in the case of non-conforming boundary conditions. Here, the modifications only concerns the relevant particle-based condition and may be implementation-specific. Hence, it is better for this to be done within the condition itself, especially since the current development of another particle Dirichlet condition based on Langrange multipliers would otherwise further complicate the code in the rotation utility.

To do so, however, requires access to the Rotate method of the utility instance from all relevant condition objects. This is not possible in the current code, since this utility class is not static (because the base class in the core is not static) and the instance is currently only available as a member of the MPMResidualBasedBossakScheme object.

This PR presents a possible way of providing access via a static member pointer in MPMParticleBaseCondition to the MPMBoundaryRotationUtility instance held by the scheme. This way, the pointer just needs to be set once in the Initialize method call of the scheme to be made available to all condition instances.

What are your thoughts on this? Are there any better ways that this can be done?

@gzjing gzjing added Discussion Proposal Suggestion of a concept for further development in a new field of application Refactor When code is moved or rewrote keeping the same behavior labels Feb 5, 2024
@VeronikaSinger
Copy link
Contributor

Thanks Guozhen for creating this PR to discuss the topic. We like to refactor/clean the boundary rotation utility in MPM application and therefore move the rotations needed for the particle based conditions locally to the conditions. This is the idea we came up with.
We are happy for any feedback or other suggestions.
@ncrescenzio , @philbucher, @sunethwarna

@VeronikaSinger
Copy link
Contributor

Hi @philbucher hi @sunethwarna
do you think the pointer to the rotation utility in the particle base condition is a suitable solution to do the rotations locally on the condition level or are there any other suggestions?

@philbucher
Copy link
Member

I dont recommend static here, since this will be shared among ALL instances of this class. Probably ok for standard sims, but perhaps one would want different ones in the future

Instead I recommend to check out how NEIGHBOR_ELEMENTS etc are handled in the fluid-app. Perhaps you can introduce a new variable that can point to the rotation utility that should be used for a particular entity. Ofc needs to be set upfront

@philbucher
Copy link
Member

Check for NEIGHBOUR_ELEMENTS

what do you think @rubenzorrilla ? Seems like a similar application to me

@gzjing
Copy link
Contributor Author

gzjing commented Mar 1, 2024

Thanks @philbucher for the feedback!

I tried out your suggestion; unfortunately, a similar implementation to NEIGHBOUR_ELEMENTS turned out to be difficult for the rotation utility because it cannot be meaningfully serialized, as its base rotation class has all of its members marked const. This prevented the definition of its pointer as an application variable, unless I define a dummy serialization for the rotation utility.

I understand your concern with the use of static, though I suppose a scenario where different rotation utilities are required would be unlikely, since that implies the DoFs and/or dimensions are different in different parts of the simulation. The only other solution I can think of would be to make the pointer an instance member, and to call a member function on the relevant Condition instance to set the pointer in Initialize(). Since this function is not present in the base class, this would necessitate the use of dynamic_cast.

What are your thoughts on this? Do you have any other suggestions?

@philbucher
Copy link
Member

I see, this is indeed a problem

Can you post the error that you get when trying to add the serialization? I was not aware that this is not working, perhaps this needs a more global change.

Other than that I also dont have good suggestions. It is indeed static vs dynamic_cast

@gzjing
Copy link
Contributor Author

gzjing commented Mar 4, 2024

The issue lies with the const Variable<double>& mrFlagVariable member of the rotation utility, which causes a -fpermissive error to be reported (namely that qualifiers would be discarded in the load method) when attempting to save/load this member.

However, even if this had worked, I'm personally not very comfortable with the idea since the current implementation of the rotation utility (with all members as const) suggests that it can neither be copied nor moved. There's two main issues I see here:

  1. Serialization requires the definition of a default constructor, which makes no sense for the rotation utility since its functionality depends on the proper specification of its const members, including a reference member that cannot be reseated later on.
  2. It feels wrong to load changes to modify a const members of the class. I realized that the serializer is using a const_cast to cast away the constness for some types, but I'm concerned that this might lead to undefined behaviour. The IMB C++ documentation suggests that "if you cast away the constness of an object that has been explicitly declared as const, and attempt to modify it, the results are undefined".

Perhaps a better solution would be to change the implementation of the rotation utility (including the base class) instead, i.e. to remove the const qualifications on the members and instead relying on ensuring all its methods are marked const.

Otherwise, I'm more inclined to the static approach for its simplicity (I'll leave a comment in the code linking to this discussion, so that the neccessary changes can be made in the future should it ever become neccessary).

Copy link
Member

@philbucher philbucher left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some general recommendations

}

// Checking whether it is normal element or penalty element
bool IsParticleBased(GeometryType& rGeometry) const
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Input should be const, otherwise you might accidentally allocate some data

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Got it!


return is_penalty;
}
bool IsParticleBased(NodeType& rNode) const
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same

{
if(this->IsSlip(rNode) )
{
const double identifier = rNode.FastGetSolutionStepValue(mrFlagVariable);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The usage of flag variables is outdated in favor of actual flags (like SLIP)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since there's no existing Kratos flag that corresponds to this use case I guess I should define a new local flag via KRATOS_DEFINE_LOCAL_APPLICATION_FLAG ?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you could do that, but flags are a bit delicate in their usage (local and global flags that use the same bit will conflict with each other)

You could still use similar code, but I recommend a variable of type bool and saved in the non-historical data. This can be accesses with GetValue. The double in your current case is a leftover from when variables could not have other types

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Understood, thanks!

}

if (CalculateStiffnessMatrixFlag == true) {
MatrixType temp_matrix = ZeroMatrix(rLeftHandSideMatrix.size1(),rLeftHandSideMatrix.size2());
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Allocating such a large matrix on the fly can have a performance impact.
Maybe you can reuse an existing variable or modify in-place?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll modify it in-place instead.

// Checking whether it is normal element or penalty element
bool IsParticleBased(GeometryType& rGeometry) const
{
for(unsigned int itNode = 0; itNode < rGeometry.PointsNumber(); ++itNode)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Better iterate over the nodes directly with geometry.Points

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Got it

@gzjing
Copy link
Contributor Author

gzjing commented Mar 25, 2024

After some discussion with @VeronikaSinger, we decided to rotate the contributions of particle-based conditions back to the global frame of reference after correcting for SLIP to avoid complications with mis-aligned frames of reference. The cleanest way to do this involves minor modifications of the base transformation utility seen in the pull request here. I hope this is okay - the changes made do not affect the existing interface and functionality as far as I'm aware.

I have also added a test case for penalty-based nonconforming slip, involving the motion of a block down a 15 degree incline with g = 10 m/s2.

I'm currently facing some issues that only crops up in Windows: there's a linking error related to the SetRotationUtility static method and its use in the scheme (see below). I've checked the relevant code, but still can't figure out where this error is coming from, especially since it doesn't show up in Ubuntu. Does anyone have an idea about how to solve this?

unity_0_cxx.obj : error LNK2019: unresolved external symbol "public: static void __cdecl Kratos::MPMParticleBaseCondition::SetRotationUtility(class Kratos::MPMBoundaryRotationUtility<class boost::numeric::ublas::matrix<double,struct boost::numeric::ublas::basic_row_major<unsigned __int64,__int64>,class boost::numeric::ublas::unbounded_array<double,class std::allocator<double> > >,class boost::numeric::ublas::vector<double,class boost::numeric::ublas::unbounded_array<double,class std::allocator<double> > > > const *)" (?SetRotationUtility@MPMParticleBaseCondition@Kratos@@SAXPEBV?$MPMBoundaryRotationUtility@V?$matrix@NU?$basic_row_major@_K_J@ublas@numeric@boost@@V?$unbounded_array@NV?$allocator@N@std@@@234@@ublas@numeric@boost@@V?$vector@NV?$unbounded_array@NV?$allocator@N@std@@@ublas@numeric@boost@@@234@@2@@Z) referenced in function "public: virtual void __cdecl Kratos::MPMResidualBasedBossakScheme<class Kratos::UblasSpace<double,class boost::numeric::ublas::compressed_matrix<double,struct boost::numeric::ublas::basic_row_major<unsigned __int64,__int64>,0,class boost::numeric::ublas::unbounded_array<unsigned __int64,class std::allocator<unsigned __int64> >,class boost::numeric::ublas::unbounded_array<double,class std::allocator<double> > >,class boost::numeric::ublas::vector<double,class boost::numeric::ublas::unbounded_array<double,class std::allocator<double> > > >,class Kratos::UblasSpace<double,class boost::numeric::ublas::matrix<double,struct boost::numeric::ublas::basic_row_major<unsigned __int64,__int64>,class boost::numeric::ublas::unbounded_array<double,class std::allocator<double> > >,class boost::numeric::ublas::vector<double,class boost::numeric::ublas::unbounded_array<double,class std::allocator<double> > > > >::Initialize(class Kratos::ModelPart &)" (?Initialize@?$MPMResidualBasedBossakScheme@V?$UblasSpace@NV?$compressed_matrix@NU?$basic_row_major@_K_J@ublas@numeric@boost@@$0A@V?$unbounded_array@_KV?$allocator@_K@std@@@234@V?$unbounded_array@NV?$allocator@N@std@@@234@@ublas@numeric@boost@@V?$vector@NV?$unbounded_array@NV?$allocator@N@std@@@ublas@numeric@boost@@@234@@Kratos@@V?$UblasSpace@NV?$matrix@NU?$basic_row_major@_K_J@ublas@numeric@boost@@V?$unbounded_array@NV?$allocator@N@std@@@234@@ublas@numeric@boost@@V?$vector@NV?$unbounded_array@NV?$allocator@N@std@@@ublas@numeric@boost@@@234@@2@@Kratos@@UEAAXAEAVModelPart@2@@Z) [D:\a\Kratos\Kratos\build\Custom\applications\MPMApplication\KratosMPMApplication.vcxproj]
D:\a\Kratos\Kratos\build\Custom\applications\MPMApplication\Custom\KratosMPMApplication.pyd : fatal error LNK1120: 1 unresolved externals [D:\a\Kratos\Kratos\build\Custom\applications\MPMApplication\KratosMPMApplication.vcxproj]

@gzjing
Copy link
Contributor Author

gzjing commented Mar 26, 2024

Thanks @philbucher :)

@gzjing gzjing marked this pull request as ready for review March 26, 2024 05:26
@gzjing gzjing requested review from a team as code owners March 26, 2024 05:26
Copy link
Member

@philbucher philbucher left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me, minor comment

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you please make a separate PR for the core changes?
This is the standard procedure, then it usually gets reviewed and merged much faster

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

PR created for core changes (#12327)

# Conflicts:
#	applications/MPMApplication/custom_strategies/schemes/mpm_residual_based_bossak_scheme.hpp
#	applications/MPMApplication/custom_utilities/mpm_boundary_rotation_utility.h
#	applications/MPMApplication/python_scripts/apply_mpm_slip_boundary_process.py
@gzjing gzjing closed this Jun 28, 2024
@gzjing gzjing reopened this Jun 28, 2024
@gzjing gzjing removed the request for review from a team July 27, 2024 06:58
@gzjing
Copy link
Contributor Author

gzjing commented Jul 27, 2024

The core changes have been merged to master in #12327 and the resulting conflicts with the changes in this PR have been resolved. This PR, which now no longer contain core changes, is now ready to be merged.

@VeronikaSinger
Copy link
Contributor

Thanks @gzjing for the changes. I tested several examples and it works as expected.

@gzjing gzjing merged commit b37159e into master Aug 9, 2024
9 of 11 checks passed
@gzjing gzjing deleted the mpm/refactor_nonconforming_slip branch August 9, 2024 18:48
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Discussion Proposal Suggestion of a concept for further development in a new field of application Refactor When code is moved or rewrote keeping the same behavior
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants