Eigensolver issues - Spurious modes & non-diagonal modal matrices

BenjaPH commented 2 years ago

I am trying to perform an eigenvalue analysis and extract all necessary modal data with the purpose of doing some external post-processes.

After seeing several examples (tests) I am confident with the eigenvalues and displacements obtained as follows,

    eigen_val_vec = model_part.ProcessInfo[KSM.EIGENVALUE_VECTOR]
    eig_vec_mat = model_part.Nodes[Id].GetValue(KSM.EIGENVECTOR_MATRIX)

but I am not sure if the following matrices contain the data I need,

    modal_mass = model_part.ProcessInfo[KSM.MODAL_MASS_MATRIX]
    modal_stiffness = model_part.ProcessInfo[KSM.MODAL_STIFFNESS_MATRIX]

I have obtained these matrices after setting in the ProjectParameters.json file the following options,

        "solver_type": "eigen_value",
        "scheme_type": "dynamic",
        "builder_and_solver_settings": {
            "use_block_builder": true
        },
        "compute_modal_decomposition": true,
        "eigensolver_diagonal_values" : {
            "mass_matrix_diagonal_value" : 1.0,
            "stiffness_matrix_diagonal_value" : 0.0
        },
        "eigensolver_settings": {
            "solver_type": "eigen_eigensystem",
            "max_iteration": 1000,
            "tolerance": 1e-06,
            "number_of_eigenvalues": 10,
            "normalize_eigenvectors": true,
            "echo_level": 1
        }

The data I expected to be in these matrices was, Modal mass ut · Mg · u Modal stiffness ut · Kg · u

As is coded in file eigensolver_strategy.hpp line 788, it seems to be that.,

    void ComputeModalDecomposition(const DenseMatrixType& rEigenvectors)
    {
        const SparseMatrixType& rMassMatrix = this->GetMassMatrix();
        SparseMatrixType m_temp = ZeroMatrix(rEigenvectors.size1(),rEigenvectors.size2());
        boost::numeric::ublas::axpy_prod(rEigenvectors,rMassMatrix,m_temp,true);
        Matrix modal_mass_matrix = ZeroMatrix(m_temp.size1(),m_temp.size1());
        boost::numeric::ublas::axpy_prod(m_temp,trans(rEigenvectors),modal_mass_matrix);
        .....
    }

Due to the orthogonality of modes both matrices should be diagonal. However, I am obtaining no null off-diagonal terms.

One thing to point out is that the model also has MPCs defined.

I would really appreciate if somebody could tell me why these matrices are not diagonal and how can I get parameters of the decoupled modal equations after projection (modal mass and modal stiffness coefficients), this is, ut · Mg · u and ut · Kg · u

Thanks in advance.

BenjaPH commented 2 years ago

Issue: In order to understand why “KSM.MODAL_MASS_MATRIX” and “KSM.MODAL_STIFFNESS_MATRIX” results sometimes in a non-diagonal matrices, we have performed a basic test with a simple supported straight beam of 10 “CrLinearBeamElement3D2N” elements.

Test: We have compared results from two identical structures.

One model with a continuous beam mesh, which is the reference model, and
the second one, with the structure broken at its mid-plane of symmetry. The two mid-span elements have been disconnected at the central point, having each side its own isolated duplicated node at centre. These two parts are rigidly joined by means of a MPC linking both coincident nodes (all d.o.f. constrained to have the same value, Ux, Uy, Uz, Rx, Ry, Rz).

A subtle detail is that the model movement is constrained to the plane OXZ.

Results: Both analyses provide the same lower eigen-frequencies and eigen-modes but only the modelization without MPCs yields diagonal “modal_mass_matrix”, “stiffness_mass_matrix” matrices. Model with MPCs gives off-diagonal terms in both matrices.

Trying to understand where is the source of this issue, I have dug into the code. May be, I have found a possible mistake, but I am a newbie in KRATOS (which is a vast quite complex program), then I would appreciate if an expert could confirm if I am right or not.

Potential mistake: With the current “ProjectParameters.json” (referred before), in the instantiation of a “class MechanicalSolver” object, the builder_and_solver selected is KratosMultiphysics.ResidualBasedBlockBuilderAndSolver.

In the source file eigensolver_strategy.hpp lines 410, 424 the master-slave constraints are applied,

this->pGetBuilderAndSolver()->ApplyConstraints(pScheme, rModelPart, rMassMatrix, b);
this->pGetBuilderAndSolver()->ApplyConstraints(pScheme, rModelPart, rStiffnessMatrix, b);

In this function, located at line 1087 of the file residualbased_block_builder_and_solver.h, it can be seen that both matrices are overwritten in the last product (mT’·A·mT),

SparseMatrixMultiplicationUtility::MatrixMultiplication(auxiliar_A_matrix, mT, rA); //A = auxilar * T   NOTE: here we are overwriting the old A matrix!

Going back to the file “eigensolver_strategy.hpp” lines 453 the “Eigenvectors” are populated with the slave dofs results

this->ReconstructSlaveSolution(Eigenvectors);

in line 788 the modal structural matrices are requested,

ComputeModalDecomposition(Eigenvectors);

This function (line 788) carry out the following operations, Eigenvectors’ · A · Eigenvectors

Remark: matrices “rMassMatrix” and “rStiffnessMatrix” were previously overwritten in the routine “ApplyConstraints” with A=(mT’·A·mT)

In my opinion, the required operation to obtain modal matrices it would be, Eigenvectors’ · A · Eigenvectors But with the original unconstrained “rMassMatrix” and “rStiffnessMatrix” matrices.

I hope someone could help to clarify this doubt.

Thanks in advance.

rubenzorrilla commented 2 years ago

Hi @BenjaPH. Thanks for the detailed explanation. Pinging the @KratosMultiphysics/structural-mechanics team as they probably know whats going on in there.

RiccardoRossi commented 2 years ago

hi @BenjaPH , i'll be even more specific and ping @armingeiser who i think is the person that knows most about this.

Now concerning the point you rise, my guess is that you are probably right, and my guess is that you are probably the first person who is looking into this specific aspect.

As you are telling there is a matrix Aorigin, but we are overwriting this matrix by Amodified = Tt Aorigin T, where the T is modified so that even the slave nodes are there as diagonal entries. The idea behind this was to avoid storing both Aorigin and Amodified for the cases in which one wanted to solve the orginal problem.

My original expectation was that

Eig' Amodified Eig

gave the same result as

Eig' Aorigin eig

but i do trust your finding that this is not the case.

I would still have a hope for the two following options

let's define an Eig_modified which is zero on all of the slave nodes and coincides with Eig_m anywhere else. Than Eig_mod' Amodified Eig_modified should hopefully give the expected result
one may also obtain the same result by doing the sum of all elemental contributions doing the Eig' A Eig element by element (with the "Eig that affects the elemental dofs)

if this it not the case than we will look into the option of not overwriting the matrix as you suggest

having said this, is there any possibility that you post the Moriginal, Koriginal and T matrices for the small problem you are dbugging with? matrix market format would be ideal

armingeiser commented 2 years ago

Hi @BenjaPH, thanks for the detailed analysis and report. I have implemented the interface to the eigen_solver, however have always just used it to obtain eigenvalues and eigenvectors of the system. As you have observed there is no issue with MPCs for these results, so i was not aware of this discrepancy for the modal matrices.

Not sure if you have tried, but since the eigenvalues and vectors are the same for both of your test cases it makes sense that with the original matrices you would also obtain the same modal matrices in the decomposition.

Out of my head i can not suggest a better solution than trying the suggestions of @RiccardoRossi and in the "worst" case store a copy of the matrices (ev. only if modal decomposition is required). It is some memory overhead of course, but compared to the memory required in the solution process it should be acceptable.

@philbucher @emiroglu FYI, i think you have worked with modal decomposition before if i am not mistaken?

Small Hint: Try if the spectra_sym_g_eigs_shift solver is faster for your use case, in general it is. https://github.com/KratosMultiphysics/Kratos/tree/master/applications/LinearSolversApplication#generalized-eigensystem-solvers

emiroglu commented 2 years ago

Hi all, @BenjaPH thanks for pointing this problem out. I remember having a discussion with @adityaghantasala that the DOFs would be rearranged during the application of MPCs. Then it is important to make sure that the matrices and the eigenvectors are matching each other during the reduction. Amodified goes with the eigenvectors before the reconstruction.

Have you had the chance to compare what is the order of magnitude of the off-diagonal values compared to the diagonal values?

Have you checked if the initial matrices are symmetric before the eigenvalue analysis? I remember that some corotational elements were resulting in unsymmetric element stiffness matrices.

I wonder how it would effect the reduction in RomApplication. Any comments from @Rbravo555 ? I remember that transient simulations were exploding with CRBeam elements + MPCs in the ROM analysis.

This would have an effect on #8174 too.

BenjaPH commented 2 years ago

Evaluation models

Please find attached the two examples that I have mentioned in my previous comment.

Continuous simple supported beam (reference model): "woMPCs_10elm.zip" woMPCs_10elm.zip
Same structure but split into two identical parts and joined with a rigid MPC: "withMPCs_10elm.zip" withMPCs_10elm.zip

You can reproduce our results just by executing the python script contained in each folder. Modal matrices are just printed out in the command windows. Projectparameters.json file with "spectra_sym_g_eigs_shift" also included.

Modal matrices

system_equations

The most natural evaluation of the modal matrices would be either with,

Complete system matrices (M, K) and expanded eigenvectors ( [T·Φ], with slaves dofs included), or
Reduced matrices (T^t·M·T, T^t·K·T) and raw eigensolution vectors (Φ, without slaves dofs)

to avoid any memory overhead i am pretty sure that the two alternatives given by @RiccardoRossi it will works fine. May be, the first one is the most straightforward option, isn’t it?.

Eigensolver results

Regarding the “quality” of the eigensolver results, there is no differences for lower modes. However, for higher ones sometimes it looses some of them, and spurious ones appear at a certain frequency. As you can see in the table below, the algorithm is not completely robust (if any, you can reproduce it with the attached models):

                        COMPARATIVE TABLE - FREQUENCY (Hz)

          |----------------- Model mesh ( 10elm ) -----------------------------|
          |------ wo MPCs ------|        |----- with MPCs ---------------------|
           (10modes)    (20modes)         (10modes)     (20modes)      (26modes)
 1         13.123056    13.123056         13.123056     13.123056      13.123056
 2         44.352929    44.352929         44.352929     44.352929      44.352929
 3         83.521319    83.521319         83.521319     83.521319      83.521319
 4         85.047338    85.047338         85.047338     85.047338      85.047338
 5        126.838223   126.838223        126.838223    126.838223     126.838223
 6        172.195671   172.195671        172.195671    172.195671     172.195671
 7        173.72397    173.723970        173.72397     173.723970     173.723970
 8        224.331591   224.331591        224.331591    224.331591     224.331591
 9        228.26061    228.260610        228.260742    228.260610     228.260610
10        248.748679   248.748679           lost       248.748679     248.748679
11                     263.572257           lost       263.572257     263.572257
12                     278.001122        278.001122    278.001122     278.001122
13                     301.248357                      301.248357     301.248357
14                     331.210482                      331.210482     331.210482
15                     361.261906                      359.136981 <*> 359.136981
16                     373.504757                      359.136981 <*> 359.136981
17                     374.477041                      359.136981 <*> 359.136981
18                     458.786549                      359.136981 <*> 359.136981
19                     466.967927                      359.136981 <*> 359.136981
20                     552.043981                      359.136981 <*> 359.136981
21                                                                    361.26190 
22                                                                    373.50475 
23                                                                    374.477041
24                                                                    458.786549
25                                                                    466.967927
26                                                                    552.043981

<*> Note that the number of spurious modes = no. of slave dofs.

Recalling the @RiccardoRossi ‘s comment

“... but we are overwriting this matrix by Amodified = Tt Aorigin T, where the T is modified so that even the slave nodes are there as diagonal entries.”

line 1491 in file residualbased_block_builder_and_solver.h

    virtual void BuildMasterSlaveConstraints(ModelPart& rModelPart)
….
        // Setting the master dofs into the T and C system
        for (auto eq_id : mMasterIds) {
            mConstantVector[eq_id] = 0.0;
            mT(eq_id, eq_id) = 1.0;
        }

        // Setting inactive slave dofs in the T and C system
        for (auto eq_id : mInactiveSlaveDofs) {
            mConstantVector[eq_id] = 0.0;
            mT(eq_id, eq_id) = 1.0;
        }

This procedure let the Mass and Stiffness matrices be of the same dimension of the global solution, with the rows and columns corresponding to the slave d.o.f. vanish (are zero), except that for the diagonal terms in both matrices (mass and stiffness) which are set on line 1116 of file residualbased_block_builder_and_solver.h to the corresponding maximum diagonals,

    void ApplyConstraints(typename TSchemeType::Pointer pScheme, ModelPart& rModelPart, TSystemMatrixType& rA, TSystemVectorType& rb) override
….
            // Apply diagonal values on slaves
            IndexPartition<std::size_t>(mSlaveIds.size()).for_each([&](std::size_t Index){
                const IndexType slave_equation_id = mSlaveIds[Index];
                if (mInactiveSlaveDofs.find(slave_equation_id) == mInactiveSlaveDofs.end()) {
                    rA(slave_equation_id, slave_equation_id) = max_diag;
                    rb[slave_equation_id] = 0.0;
                }

This is an acceptable artifice for static analyses, but I guess that in the modal eigenproblem (see equation above) introduce spurious modes of frequency = sqrt(kmax_diag/mmax_diag). In our example are just 6 modes (= 6 slave dofs) and come out at 359.136981 Hz.

Something similar is done with the DirichletConditions, restrained dofs are also preserved in the eigensolution, and their diagonal terms are scaled by a given function argument “Factor”, line 674 in file eigensolver_strategy.hpp

    void ApplyDirichletConditions(SparseMatrixType& rA, double Factor)
                // row dof is fixed. zero off-diagonal columns and factor diagonal
                for (std::size_t j = ColBegin; j < ColEnd; ++j)
                {
                    if (static_cast<int>(AColIndices[j]) != k)
                    {
                        AValues[j] = 0.0;
                    }
                    else
                    {
                        AValues[j] *= Factor;
                    }
                }

in this case, Factor=mStiffnessMatrixDiagonalValue = 1.0 and Factor=mMassMatrixDiagonalValue=0, this would result in the injection of eigenvalues at infinite. The “eigen” library algorithm it is supposed to cope with singular mass matrices, but it is for sure an additional difficulty. I have not any explanation for the modes lost, but the mass matrix singularity could be a possible cause.

By the way, mT is computed twice, one time for each call to the function ‘ApplyConstraints’. This function is called to obtain both the constrained rStiffnessMatrix and constrained rMassMatrix, 2 x BuildMasterSlaveConstraints(rModelPart);

This is my personal opinion. Please take these comments with caution. I hope this can be of any help.

PS @RiccardoRossi I am trying to write the matrices in matrix market format. It take me a while, I need to deal with the compilation …. @emiroglu I hope i can answer your questions about symmetry soon.

loumalouomega commented 2 years ago

Question, using elimination builder and solver with constraints may solve the issue?

loumalouomega commented 2 years ago

Question, using elimination builder and solver with constraints may solve the issue?

(at least the sporious part)

rubenzorrilla commented 2 years ago

Question, using elimination builder and solver with constraints may solve the issue?

As far as I know the elimination B&S is not MPC-compatible yet (maybe I'm wrong).

loumalouomega commented 2 years ago

Question, using elimination builder and solver with constraints may solve the issue?

As far as I know the elimination B&S is not MPC-compatible yet (maybe I'm wrong).

Yes, it was me who implemented (there a little OpenMp bug I cannot find), it is independent of the standard one

BenjaPH commented 2 years ago

Following your suggestions, i have done some tests with the Beam model including MPCs.

ResidualBasedBlockBuilderAndSolver with spectra_sym library: Catch the first natural frequency but the rest are spurious.
ResidualBasedEliminationBuilderAndSolverWithConstraints: give similar results to the reference model up to the 7th mode, but after this one, some spurious modes appear.
ResidualBasedBlockBuilderAndSolverWithLagrangeMultiplier: the firsts eigenfactors are negative (is not possible to obtain real frequencies), the following ones are nearly zero.


                        COMPARATIVE TABLE - FREQUENCY (Hz)         <|>   negative lambda factor
                                                                         no real sqrt()
          |------------------------------ Model mesh ( 10elm ) --------------------------|
          | wo MPCs |    |------------------------ with MPCs ----------------------------|
          (Reference)    (Rb_BB&S_eig) (Rb_BB&S_spt)  (Rb_EB&S_eig)     (Rb_BB&S_LM_eig)       
 1          13.123056       13.123056    13.123056     13.123056      -75117.91820204783  *
 2          44.352929       44.352929    13.929155*     44.352929      -21448.755365459205*
 3          83.521319       83.521319    13.929210*     83.521319   -2.161535571554726e-03*
 4          85.047338       85.047338    13.929210*     85.047338   -6.275719300479650e-04*
 5         126.838223      126.838223    13.929210*    126.838223   -4.918551342093946e-05*
 6         172.195671      172.195671    13.929210*    172.195671   -3.888269205286084e-06*
 7         173.723970      173.723970    13.929210*    173.723970   -1.255736757037714e-07*
 8         224.331591      224.331591    13.929210*    198.431527*  -1.069627659847707e-08*
 9         228.260610      228.260610    13.929210*    224.331591   -3.436626353232495e-09*
10         248.748679      248.748679    13.929210*    228.260610   -3.032532192128545e-10*
11         263.572257      263.572257    13.929210*    248.748679   -6.379973002282229e-11*
12         278.001122      278.001122    13.929210*    263.572257   -5.171982882275921e-11*
13         301.248357      301.248357    13.929210*    278.001122   -6.646883191289502e-12*
14         331.210482      331.210482    13.929210*    292.451060*  -7.441249339044151e-44*
15         361.261906      359.136981*   13.929210*    301.248357   1.1456279024815887e-49*
16         373.504757      359.136981*   13.929210*    331.210482   1.3723424980234234e-33*
17         374.477041      359.136981*   13.929210*    361.261906    5.886084936787898e-26*
18         458.786549      359.136981*   13.929210*    373.504757   1.4811905545020017e-09*
19         466.967927      359.136981*   13.929210*    374.477041   1.6980771438722432e-09*
20         552.043981      359.136981*   13.929210*    458.786549   2.5915256804442248e-08*

<*>  spurious modes

Rb_BB&S_eig -> ResidualBasedBlockBuilderAndSolver_eigen
Rb_BB&S_spt -> ResidualBasedBlockBuilderAndSolver_spectra_sym
Rb_EB&S_eig -> ResidualBasedEliminationBuilderAndSolverWithConstraints_eigen
Rb_BB&S_LM_eig -> ResidualBasedBlockBuilderAndSolverWithLagrangeMultiplier_eigen ()

BenjaPH commented 2 years ago

Just to close this issue from my side, please find attached the system matrices for the testing problem with MPCs sent on the 7th of April (requested by @RiccardoRossi ). As can be seen, global matrices M, K are fully symmetric (answering the question of @emiroglu).

System-Matrices.zip

Working out the system matrices - following the conventional approach - with a kinematic relation matrix 'T' relating the complete set of dofs with only the master + free dofs (rectangular matrix instead of that used in Kratos, squared matrix).

and subsequently removing rows and columns corresponding to restrained dofs, the eigenvalues fully agree with the example without MPCs.

Of course, the conventional approach requires more memory to hold simultaneously the reduced matrices but yield correct results.

It seems that the method implemented in Kratos introduce spurious eigenvalues when the model contains MPCs.

RiccardoRossi commented 2 years ago

sorry, i would prefer to keep it open. I have this in mind and eventually i will try to have this solved "well"

BenjaPH commented 2 years ago

Great !!!

loumalouomega commented 2 years ago

If there is a bug in the MPC I am part responsible @RiccardoRossi

BTW: https://scicomp.stackexchange.com/questions/31260/solver-for-generalized-eigenvalue-problem-with-multipoint-constraints

BenjaPH commented 2 years ago

Is there any update on this issue?

I would be very glad to know about any progress …

loumalouomega commented 2 years ago

Is there any update on this issue?

I would be very glad to know about any progress …

The problem was MPC + EigenValues?¿

loumalouomega commented 2 years ago

Is there any update on this issue? I would be very glad to know about any progress …

The problem was MPC + EigenValues?¿

Okay, after reading the previous messages, I see that the problem is that we don't apply MPC to the mass matrix assembled...

philbucher commented 2 years ago

Is there any update on this issue? I would be very glad to know about any progress …

The problem was MPC + EigenValues?¿

Okay, after reading the previous messages, I see that the problem is that we don't apply MPC to the mass matrix assembled...

we do: https://github.com/KratosMultiphysics/Kratos/blob/6c1557eddc1d86f2e7f02d3b4a7efffeef568ece/applications/StructuralMechanicsApplication/custom_strategies/custom_strategies/eigensolver_strategy.hpp#L409-L411

loumalouomega commented 2 years ago

Is there any update on this issue? I would be very glad to know about any progress …

The problem was MPC + EigenValues?¿

Okay, after reading the previous messages, I see that the problem is that we don't apply MPC to the mass matrix assembled...

we do:

https://github.com/KratosMultiphysics/Kratos/blob/6c1557eddc1d86f2e7f02d3b4a7efffeef568ece/applications/StructuralMechanicsApplication/custom_strategies/custom_strategies/eigensolver_strategy.hpp#L409-L411

In that case I dont know what to do

KratosMultiphysics / Kratos