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

[Core] Neighbor MPC Assigner: A process and utility for assigning MPCs to neighbouring nodes #11083

Merged
merged 28 commits into from
May 9, 2023

Conversation

SADPR
Copy link
Contributor

@SADPR SADPR commented May 5, 2023

DESCRIPTION:

Utility:

This pull request adds a new class that provides a method to search for neighboring nodes of a given nodes structure within a specified radius and assign a multipoint constraint (MPC). The MPC weights are calculated based on a spatial function that employs radial basis functions.

Process:

This pull request enhances the functionality of the master-slave MPC algorithm by facilitating the discovery of neighboring nodes in a master model part within a designated radius for each node in the slave model part. The process involves the following steps:

  1. For each node in the slave model part, the utility searches for neighboring nodes in the master model part within a specified radius.
    image
    image
    image

  2. The utility then establishes an MPC between the slave node and its neighboring master nodes, and calculates the MPC weights based on a spatial function that employs radial basis functions.
    image
    image

  3. The process repeats for all nodes in the slave model part.

Added files:

  • assign_mpcs_to_neighbours_utility.cpp: This file contains all the functions to look for neighboring nodes and apply multipoint constraints (MPCs) using radial basis functions.
  • assign_mpcs_to_neighbours_utility.h: This header file defines the functions and data structures used in assign_mpcs_to_neighbours_utility.cpp`.
  • assign_mpcs_to_neighbours_process.py: This Python script calls the functions in assign_mpcs_to_neighbours_utility.cpp in the appropriate order to assign MPCs to neighboring nodes.
  • assign_mpcs_to_neighbours_process.mdpa: This is an MDPA file used for testing the MPC algorithm.

Modified files:

  • add_geometrical_utilities_to_python.cpp: This file was modified to export the functions in assign_mpcs_to_neighbours_utility.cpp to Python, so that they can be called from assign_mpcs_to_neighbours_process.py.
  • test_processes.py: This file was modified to include a small test to validate the functionality of assign_mpcs_to_neighbours_process.py.

SADPR added 4 commits May 3, 2023 11:32
@SADPR SADPR requested a review from a team as a code owner May 5, 2023 09:01
@loumalouomega
Copy link
Member

@KratosMultiphysics/altair and specially @ddiezrod

@loumalouomega loumalouomega changed the title Neighbor MPC Assigner: A process and utility for assigning MPCs to neighboring nodes [Core] Neighbor MPC Assigner: A process and utility for assigning MPCs to neighbouring nodes May 5, 2023
Comment on lines 169 to 170
for constraint in self.computing_model_part.MasterSlaveConstraints:
constraint.Set(KM.TO_ERASE)
Copy link
Member

Choose a reason for hiding this comment

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

Use VariableUtils

@@ -0,0 +1,3308 @@
Begin Properties 0
Copy link
Member

Choose a reason for hiding this comment

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

Please, can you consider a smaller mesh?, you can put this one in the Example repository


ModelPart::MasterSlaveConstraintType const& r_clone_constraint = KratosComponents<MasterSlaveConstraint>::Get("LinearMasterSlaveConstraint");

#pragma omp parallel
Copy link
Member

Choose a reason for hiding this comment

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

Please ParallelUtilities

DofPointerVectorType& rCloud_of_dofs,
array_1d<double,3>& rSlave_coordinates
){
KRATOS_TRY;
Copy link
Member

@loumalouomega loumalouomega May 5, 2023

Choose a reason for hiding this comment

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

Tabs missing (spaces)

Matrix& rCloud_of_nodes_coordinates
)
{
KRATOS_TRY;
Copy link
Member

Choose a reason for hiding this comment

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

Tabs missing (spaces)

)
{
//TODO: Do it in parallel
KRATOS_TRY;
Copy link
Member

Choose a reason for hiding this comment

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

Tabs missing (spaces)

Copy link
Member

@loumalouomega loumalouomega left a comment

Choose a reason for hiding this comment

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

Some comments

@SADPR
Copy link
Contributor Author

SADPR commented May 5, 2023

I will use a smaller mesh for the test and update to Parallel Utilities.
I will upload this case to the Kratos Examples.
Assign_MPCs_to_neighbours

@loumalouomega
Copy link
Member

I cannot unsee:
image

image

Matrix r_cloud_of_nodes_coordinates;
array_1d<double,3> r_slave_coordinates;
GetDofsAndCoordinatesForNode(pNode, rVariable, r_slave_dof, r_slave_coordinates);
KRATOS_WATCH(i)
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
KRATOS_WATCH(i)

Debugging leftover.

@SADPR
Copy link
Contributor Author

SADPR commented May 7, 2023

I've been trying to switch from OMP to block_for_each the implementation, and I did the following tests to have a counter "i" that could be thread-safe. The reason for using the "i" counter is to avoid repeating an MPC's Id.

  • Using std::atomic<int>: This test showed repeating values in some threads, which means it's not thread-safe in this case.
  • Using std::atomic_flag: This test showed repeating values in some threads as well, indicating it's not thread-safe either.
  • Using a std::mutex with std::lock_guard: This test worked correctly but was around 50 times slower than using pragma omp.
  • Using a std::mutex with std::unique_lock: This test did not work, as it resulted in repeating values in some threads.
  • Using #pragma omp parallel for: This test worked correctly and was faster than using std::mutex with std::lock_guard.

Overall, it seems like using std::mutex with std::lock_guard is the only method that worked correctly for this particular use case, but it comes with a significant performance penalty. The pragma omp approach was faster, but you mentioned you wanted to avoid using OpenMP.

For the current implementation, I kept the original code using #pragma omp in the function AssignMPCsToNodes, and I added the implementation with block_for_each (std::mutex and std::lock_guard) in the function AssignMPCsToNodesParallelUtilities. Additionally, I have tried to use the KRATOS_CRITICAL_SECTION as a replacement for the #pragma omp critical section in order to add the MPCs thread-safely, but I haven't been successful yet. I am still looking for a way to substitute the OMP critical section for adding the MPCs thread-safely. If you can point me in the right direction, it would be greatly appreciated.

* @param RotationMatrix 3x3 rotation matrix to be applied to the nodes' positions.
**/
void AssignRotationToNodes(
NodesContainerType pNodes,
Copy link
Member

Choose a reason for hiding this comment

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

This should be a reference. Indeed it is strange to me that this works if you pass it by copy...

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Not used anymore.


// Search for the neighbouring nodes (in rStructureNodes) of each rNode on a given array of rNodes
void SearchNodesInRadiusForNodes(
NodesContainerType const& rNodes,
Copy link
Member

Choose a reason for hiding this comment

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

Though it is valid, in Kratos we normally put the const in front.


// Search for the neighbouring nodes (in rStructureNodes) of each rNode on a given array of rNodes
void AssignMPCsToNeighboursUtility::SearchNodesInRadiusForNodes(
NodesContainerType const& rNodes,
Copy link
Member

Choose a reason for hiding this comment

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

Normally in Kratos we put the const in front.

Comment on lines +53 to +54
VectorResultNodesContainerType& rResults,
const double MinNumOfNeighNodes)
Copy link
Member

Choose a reason for hiding this comment

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

It is a good programming practice to put the results arguments after the input (const ones).

const double MinNumOfNeighNodes)
{
KRATOS_TRY;
NodesContainerType::ContainerType& nodes_array = const_cast<NodesContainerType::ContainerType&>(rNodes.GetContainer());
Copy link
Member

Choose a reason for hiding this comment

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

auto& r_nodes_array

Copy link
Member

Choose a reason for hiding this comment

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

And why you need the const_cast?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Resolved, changed to:

    KRATOS_TRY;
    auto& r_nodes_array = rNodes.GetContainer();

    rResults.resize(r_nodes_array.size());

    // Declare a counter variable outside the loop as std::atomic<int>
    std::atomic<int> i(0);

    block_for_each(r_nodes_array, [&](const Node<3>::Pointer& pNode)
    {


block_for_each(nodes_array, [&](Node<3>::Pointer& pNode)
{
double localRadius = Radius;
Copy link
Member

Choose a reason for hiding this comment

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

Please follow the style guide (local_radius).

KRATOS_TRY;
NodesContainerType::ContainerType& nodes_array = const_cast<NodesContainerType::ContainerType&>(rNodes.GetContainer());

rResults.resize(nodes_array.size());
Copy link
Member

Choose a reason for hiding this comment

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

This means you're allocating much memory than needed. I'd do a shrink_to_fit at the end of this function to avoid so.

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 am allocating what I need. Nonetheless, I could updated it to this:

    // Resize rResults vector
    rResults.resize(r_nodes_array.size());

    // Declare a counter variable outside the loop as std::atomic<int>
    std::atomic<int> i(0);

    block_for_each(r_nodes_array, [&](const Node<3>::Pointer& pNode)
    {
        double local_radius = Radius;

        IndexType it = i.fetch_add(1); // Atomically increment the counter and get the previous value

        while (rResults[it].size()<MinNumOfNeighNodes){
            rResults[it].clear();
            SearchNodesInRadiusForNode(pNode, local_radius, rResults[it]);
            local_radius += Radius;
        }
    });

    // Optimize memory usage
    rResults.shrink_to_fit();

Copy link
Member

Choose a reason for hiding this comment

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

No. You are allocating much more than what you need since you're resizing the vector to the number of nodes, which would only be the case if all the nodes are found. Indeed, the proper usage would be

rResults.reserve(r_nodes_array.size());

// ... your stuff in here ...

rResults.shrink_to_fit();

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 am changing the names of the functions to make them more self-explanatory. In this function, we retrieve the degrees of freedom (DOFs) and coordinates of the cloud of nodes. Therefore, we already know the size.


//Search for the neighbouring nodes (in rStructureNodes) of the given rNode
void AssignMPCsToNeighboursUtility::SearchNodesInRadiusForNode(
NodeType::Pointer pNode,
Copy link
Member

Choose a reason for hiding this comment

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

I'd pass a constant reference instead. Note that the SearchInRadius supports it.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Do you find it necesarry, since this is a constant pointer. Isn't it okay like this?

Copy link
Member

Choose a reason for hiding this comment

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

I don't think it is a constant pointer but nevermind, leave it like this.


if (ComponentIndex < 3) {
std::string component_variable_name = rVectorVariable.Name() + component_suffixes[ComponentIndex];
return KratosComponents<Variable<double>>::Get(component_variable_name);
Copy link
Member

Choose a reason for hiding this comment

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

This is extremely expensive and being done, if I'm not wrong, inside a loop. An alternative approach must be used.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

After a long afternoon, I managed to do it. Now it only acces the Kratos components once and has a pointer in the loop.

Comment on lines +121 to +122
rCloud_of_dofs.push_back(pNode->pGetDof(rVariable));
rSlave_coordinates = pNode->Coordinates();
Copy link
Member

Choose a reason for hiding this comment

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

I don't see the point on doing a function of two lines. IMO it only makes more complex browsing the code with no advantage at all.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It is to preserve similarity. Is it really necessary to undo this function? I know it is not typical to have a two line function, but for me, it preserves similarity.

Copy link
Member

Choose a reason for hiding this comment

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

I don't get your point on "similarity"... I'd remove it but if you really require it leave it. Up to you.

Comment on lines +138 to +140
if (pNode->HasDofFor(GetComponentVariable(rVariable, 0)) &&
pNode->HasDofFor(GetComponentVariable(rVariable, 1)) &&
pNode->HasDofFor(GetComponentVariable(rVariable, 2))) {
Copy link
Member

Choose a reason for hiding this comment

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

What would happen in the 2D case in which the third component DOF might not be present?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Is this really something that happens?
The user can also specify ["DISPLACEMENT_X","DISPLACEMENT_Y"].

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 thought the 2D also 3 components. But solves only for X and Y components.

Copy link
Member

Choose a reason for hiding this comment

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

Depends on the implementation so you cannot take it as granted.

DofPointerVectorType& rCloud_of_dofs_x,
DofPointerVectorType& rCloud_of_dofs_y,
DofPointerVectorType& rCloud_of_dofs_z,
array_1d<double,3>& rSlave_coordinates
Copy link
Member

Choose a reason for hiding this comment

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

Please follow the style guide for function arguments.



void AssignMPCsToNeighboursUtility::GetDofsAndCoordinatesForNodes(
ResultNodesContainerType nodes_array,
Copy link
Member

Choose a reason for hiding this comment

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

Why passing by copy?

Comment on lines +182 to +185
rCloud_of_dofs_x.resize(nodes_array.size());
rCloud_of_dofs_y.resize(nodes_array.size());
rCloud_of_dofs_z.resize(nodes_array.size());
rCloud_of_nodes_coordinates.resize(nodes_array.size(), 3);
Copy link
Member

Choose a reason for hiding this comment

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

I'd try to shrink these containers afterwards.

Copy link
Member

Choose a reason for hiding this comment

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

(same reseve and shrink_to_fit comment).

// Check if the node has the required DOFs
if (pNode->HasDofFor(GetComponentVariable(rVariable, 0)) &&
pNode->HasDofFor(GetComponentVariable(rVariable, 1)) &&
pNode->HasDofFor(GetComponentVariable(rVariable, 2))) {
Copy link
Member

Choose a reason for hiding this comment

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

Same comment about the 2D third component.

{
//TODO: Do it in parallel
KRATOS_TRY;
NodesContainerType::ContainerType& nodes_array = const_cast<NodesContainerType::ContainerType&>(pNodes.GetContainer());
Copy link
Member

Choose a reason for hiding this comment

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

auto& r_nodes_array

Copy link
Member

Choose a reason for hiding this comment

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

Why the const_cast?


// Assign rotation to given nodes
for(unsigned int i = 0; i < nodes_array.size(); i++){
auto coordinate_vector = nodes_array[i]->GetInitialPosition().Coordinates();
Copy link
Member

Choose a reason for hiding this comment

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

You are allocating/deallocating an unnecessary copy for each node.

// Assign rotation to given nodes
for(unsigned int i = 0; i < nodes_array.size(); i++){
auto coordinate_vector = nodes_array[i]->GetInitialPosition().Coordinates();
Vector rotated_position(coordinate_vector.size());
Copy link
Member

Choose a reason for hiding this comment

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

Can't this be an array_1d<double,3>?

Comment on lines +222 to +224
nodes_array[i]->X() = rotated_position[0];
nodes_array[i]->Y() = rotated_position[1];
nodes_array[i]->Z() = rotated_position[2];
Copy link
Member

Choose a reason for hiding this comment

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

noalias(nodes_array[i]->Coordinates()) = rotated_position; I think this should work.

Comment on lines +302 to +308
void AssignMPCsToNeighboursUtility::AssignMPCsToNodes(
NodesContainerType pNodes,
double const Radius,
ModelPart& rComputingModelPart,
const Variable<array_1d<double, 3>>& rVariable,
double const MinNumOfNeighNodes
)
Copy link
Member

Choose a reason for hiding this comment

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

This is basically copy & paste of the above function but modifying the section in which the Master-Slave Constraints are created. I'd do the creation within a private auxiliary function and have a unique implementation of AssignMPCsToNodes.

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 have a concern here, at first I though of the user passing a variable which could be either array_1d or double.
If array_1d, then you have to create the variable components lists, etc. If double, you can directly retrieve the varible with GetDof. Now I am thinking to always create a list of variables, e.g, if ["VELOCITY","PRESSURE], then create a list:
["DISPLACEMENT_X","DISPLACEMENT_Y","PRESSURE"] for 2D.
Or:
["DISPLACEMENT_X","DISPLACEMENT_Y","DISPLACEMENT_Z","PRESSURE"] for 3D.
We can discuss this in person, since it is taking a lot of time.

Copy link
Member

Choose a reason for hiding this comment

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

OK.

// Add the constraints to the rComputingModelPart in a single call
rComputingModelPart.AddMasterSlaveConstraints(temp_constraints.begin(), temp_constraints.end());

KRATOS_INFO("AssignMPCsToNeighboursUtility") << "Build and Assign MPCs Time: " << build_and_assign_mpcs.ElapsedSeconds() << std::endl;
Copy link
Member

Choose a reason for hiding this comment

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

This may result in a too verbose output. If you want to output this, I'd introduce an echo level.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

TODO.


BuiltinTimer build_and_assign_mpcs;

int prev_num_mpcs = rComputingModelPart.NumberOfMasterSlaveConstraints();
Copy link
Member

Choose a reason for hiding this comment

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

In here you are assuming that your new MPCs are not interacting with previous ones. Though this might be a valid assumption, I'd throw a warning message if the previous number of MPCs is non-zero to inform about this.


int prev_num_mpcs = rComputingModelPart.NumberOfMasterSlaveConstraints();

NodesContainerType::ContainerType& nodes_array = const_cast<NodesContainerType::ContainerType&>(pNodes.GetContainer());
Copy link
Member

Choose a reason for hiding this comment

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

auto& r_nodes_array.


NodesContainerType::ContainerType& nodes_array = const_cast<NodesContainerType::ContainerType&>(pNodes.GetContainer());

ModelPart::MasterSlaveConstraintType const& r_clone_constraint = KratosComponents<MasterSlaveConstraint>::Get("LinearMasterSlaveConstraint");
Copy link
Member

Choose a reason for hiding this comment

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

const auto& r_clone_constraint`.

Comment on lines +268 to +270
DofPointerVectorType r_cloud_of_dofs, r_slave_dof;
Matrix r_cloud_of_nodes_coordinates;
array_1d<double, 3> r_slave_coordinates;
Copy link
Member

@rubenzorrilla rubenzorrilla May 10, 2023

Choose a reason for hiding this comment

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

You're allocating/deallocating these for each node. Declare them as a TLS outside the loop.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Master-Slave constraints require rSlaveDofs to be a DofPointerVectorType. To do r_slave_dofs thread_local I need to resize and assign to Dof::Pointer. Or:
rCloudOfDofs.resize(1);
rCloudOfDofs[0] = pNode->pGetDof(rVariable);
Which I dont like. I will fix it.

Comment on lines +2559 to +2562
list_of_processes = process_factory.KratosProcessFactory(current_model).ConstructListOfProcesses(settings["process_list"])
for process in list_of_processes:
process.ExecuteInitialize()
process.ExecuteInitializeSolutionStep()
Copy link
Member

Choose a reason for hiding this comment

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

You have a unique process. I don't see the point on doing it with a list. I'd directly call the constructor.

Comment on lines +2593 to +2594
for ew, ow in zip(expected_weights, obtained_weights):
self.assertAlmostEqual(ew, ow, delta=1e-6) # Assert the weights are equal within the specified tolerance
Copy link
Member

Choose a reason for hiding this comment

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

If I'm not wrong, there is an assertAlmostEqualVector or something of this sort.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

assertVectorAlmostEqual

Copy link
Member

@rubenzorrilla rubenzorrilla left a comment

Choose a reason for hiding this comment

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

Though this is a working implementation, I think it doesn't suffice the quality requirements for being in the KratosCore. I hope you take my comments into consideration in a new PR.

KRATOS_CATCH("");
}

void AssignMPCsToNeighboursUtility::AssignRotationToNodes(
Copy link
Member

@rubenzorrilla rubenzorrilla May 10, 2023

Choose a reason for hiding this comment

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

This is not used anywhere in this process. Besides, it is a custom function for your target application. Long story short, it must be removed from here to be placed somewhere else.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants