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] Incorrect AMGCLSolver::mBlockSize #11853

Open
matekelemen opened this issue Nov 25, 2023 · 5 comments
Open

[Core] Incorrect AMGCLSolver::mBlockSize #11853

matekelemen opened this issue Nov 25, 2023 · 5 comments
Assignees

Comments

@matekelemen
Copy link
Contributor

matekelemen commented Nov 25, 2023

Judging by how it's calculated, I'm assuming AMGCLSolver::mBlockSize should represent the number of DoFs per variable (or node, I'm not sure), but it seems to me that the current implementation is broken.

Can someone explain what exactly is mBlockSize supposed to be equal to on the Kratos side?

I'm running a linear structural mechanics case with quadratic elements and some MPCs. For this kind of model, the number of DoFs in a variable should be 3 (x, y, z displacement or rotation components). However, AMGCLSolver ends up with an mBlockSize of 1, which completely derails the solver to the point it doesn't even converge (hardcoding it to 3 results in convergence).

Here's how mBlockSize is calculated:

unsigned int old_node_id = rDofSet.size() ? rDofSet.begin()->Id() : 0;
for (auto it = rDofSet.begin(); it!=rDofSet.end(); it++) {
if(it->EquationId() < TSparseSpaceType::Size1(rA) ) {
IndexType id = it->Id();
if(id != old_node_id) {
old_node_id = id;
if(old_ndof == -1) old_ndof = ndof;
else if(old_ndof != ndof) { //if it is different than the block size is 1
old_ndof = -1;
break;
}
ndof=1;
} else {
ndof++;
}
}
}
if(old_ndof == -1)
mBlockSize = 1;
else
mBlockSize = ndof;

It appears to me that this piece of code is supposed to count how many DoFs share the same ID, and makes 2 assumptions:

  1. every ID is referenced by the same number of DoFs, so it's enough to count once (the first ID)
  2. DoFs sharing their IDs are consecutive in rDofSet

The first assumption seems reasonable to me, but the second does not. rDofSet is a PointerVectorSet<Dof<double>>, which means that it's sorted by the values of raw Dof pointers. I don't see how this would guarantee consecutive IDs.

Am I missing something or is this a bug?

P.S.: the solver does actually converge with a block size of 1, but more than 500% slower (for my specific model).

@AlejandroCornejo
Copy link
Member

Wow, that's very interesting (as much as concerning...). I had the same idea about what is the Block size and seems that this piece of code is not workling as expected... I pinned @ddemidov and @RiccardoRossi to be sure about this

@ddemidov
Copy link
Member

I am not familiar with Kratos internals well enough to understand what is the problem here, sorry.

@RiccardoRossi
Copy link
Member

RiccardoRossi commented Nov 29, 2023

for (aggregation) AMG to work correctly you need blocks of dofs to be repeated for every node.
let's imagine that we have dx dy dz to be the displacements

if we have

dx_1 dy_1 dz_1 dx_2 dy_2 dz_2 dx_3 dy_3 dz_3

than all of the nodes 1 2 and 3 have the same anount of dofs. 3 in this case.

however if you have for example

dx_1 dy_1 dz_1 dx_2 dy_2 dx_3 dy_3 dz_3

that is, for example if 2 does not have "dz" in the system, BlockSize will be detected as 1 as it is the only multiplier we have

the "BlockBuilderAndSolver" ensures that this property is kept. The Eliminiation builder and solver does not ...

to be clear, if you have a multidimensional problem (like structures) and the BlockSize==1 aggregation multigrid will NOT work correctly. You could still use RugenStuben for example, but aggregation based amg will not work properly

@dcagritan
Copy link
Contributor

Do we have documentation regarding solvers in Kratos? Documentation or a wiki page that explains the parameters of the solvers for instance?

@RiccardoRossi
Copy link
Member

RiccardoRossi commented Dec 2, 2023 via email

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

No branches or pull requests

6 participants