Stiffness and mass matrix in Python

[MagaFem]

Is possible to get the mass- and stiffness-matrix to python-workspace for a structure mechanical model? Looks like matrices are assembled in: eigensolver_strategy.hpp

Thanks!

[loumalouomega]

Hi, there are public methods in C++, we can add it easily to python, but right now is not a good moment, we are migrating our python interface. You can change your code locally

[loumalouomega]

Go to add_custom_strategies_to_python in StructuralMechanicsApplication and add change the EigenSolver declaration as follows:

    // Eigensolver Strategy
    class_< EigensolverStrategyType, bases< BaseSolvingStrategyType >, boost::noncopyable >
            (
                "EigensolverStrategy", init<ModelPart&, BaseSchemeType::Pointer, BuilderAndSolverPointer>() )
                .def("GetMassMatrix", &EigensolverStrategyType::GetMassMatrix)
                .def("GetStiffnessMatrix", &EigensolverStrategyType::GetStiffnessMatrix)
            ;

[RiccardoRossi]

@qaumann or @armingeiser is this solution fine? shall a PR be done to add this capability?

[armingeiser]

I think this could be useful, if one wants to solve the eigenvalue problem externally or using e.g. numpy. But one should keep in mind that the Eigensolver Strategy does modify the main diagonal entries of the supported dofs.

Edit: using Eigens capabilities, one could even translate them to numpy matices easily (as far as i know)

[qaumann]

@MagaFem, what are you planning to do with the system matrices?

For me it seems @MagaFem is interested in the mass or stiffness matrices in general, not only for eigenvalue calculations. So it could be nice to add an "export to python" functionality for the system matrices, not limited to the eigenvalue case and thus export the unmodified matrices ignoring the boundary conditions.

[RiccardoRossi]

@qaumann i agree that capability is interesting, however the problem is that the only one that could do this is the builder and solver.

One could achieve the desired behaviour by calling the builder and solver directly. One important problem is that it is also needed to get the doflist, to be able to know how ids are associated to nodes and dofs.

[RiccardoRossi]

we'll propose an interface to do this in a useful way

[RiccardoRossi]

ok, so this is the first iteration:

template< TSparseSpace >
class MatrixBuilder
{
      void MatrixBuilder(rModelPart, Flag BuilderAndSolverType=BlockBuilderAndsolver)
      void Initialize()
      void Clear()

      #NOTE: one may change behaviour by passing variables through rModelPart.ProcessInfo if needed
      TSparseSpace::VectorType&  ComputeRHS( bool ApplyDirichletConditions )
      TSparseSpace::SparseMatrixType& ComputeStiffnessMatrix( bool ApplyDirichletConditions )
      TSparseSpace::SparseMatrixType& ComputeMassMatrix( bool ApplyDirichletConditions )
      TSparseSpace::SparseMatrixType& ComputeDampingMatrix( bool ApplyDirichletConditions )

      ModelPart::DofsArrayType& GetDofList()

}

we shall also complete this by an interface to construct scipy matrices from our matrices, but this is easily done on the python side.

Also we should improve the interface of the ModelPart::DofsArrayType so to be able to:

1. iterate over dofs in python
2. get the variable associated
3. get if fixed or not

[philbucher]

+1 for the idea

One thought since you want to construct scipy matrices: I think pybind can also work with Eigen directly, so im am putting @armingeiser and @oberbichler in the loop since they might be interested

[oberbichler]

PyBind has build-in converters for Numpy/SciPy <-> Eigen. The exchange can be done without copying the underlaying data. It also provides a lot of helpers for inherit C++ classes in Python.

During my masterthesis I used this interface e.g. for developing iga-elements in python (maximum comfort) and put them in the C++ model (maximum speed). It works very well out-of-box.

I don't knot how difficult it is to write equivalent wrappers for ublas/Kratos but I think they would be very usefull!

[RiccardoRossi]

Hi @oberbichler i think that we can quite easily export the csr arrays to numpy using pybind. once that is done scipy has a constructor with that.

@philbucher , I would prefer not to depend on eigen for this (recall that Eigen needs to be strictly optional)

[MagaFem]

@qaumann right, I had a couple things in mind, namely model reduction and calculation of compliance. I saw that the eigensolver_strategy.hpp assembled and named the matrices and that is why I referred to the eigenvalue case.

[RiccardoRossi]

Is there anyone volunteering to implement this? (or the changes we agree)

[loumalouomega]

Do you mean the interface I suggested at the beginning or something more complete?

[loumalouomega]

The interface I suggested will export to Kratos Vector and Matrix, for numpy some more extra work:

https://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html

[RiccardoRossi]

The MatrixBuilder i am suggesting in one of the lates posts

[loumalouomega]

"Hi @oberbichler i think that we can quite easily export the csr arrays to numpy using pybind. once that is done scipy has a constructor with that. " This one I assume, do we add as an independent utility or do you prefer to add directly in add_matrix_to_python on the core?

[RiccardoRossi]

closing as inactive. Reopen if needed or when someone wants to volunteer for the implementation

[loumalouomega]

@RiccardoRossi do you mean this? https://pybind11.readthedocs.io/en/stable/advanced/cast/eigen.html

[veiguf]

This seems to be solved, but I do not understand how I can get the stiffness and mass matrix to Python. Is it documented anywhere? By the way: Is there an easy and clean way to convert a Kratos.Matrix to a Numpy array?

[loumalouomega]

This seems to be solved, but I do not understand how I can get the stiffness and mass matrix to Python. Is it documented anywhere?

The LHS is easy, the solvers could be adapted in order to retrieve the K and M in separated way

By the way: Is there an easy and clean way to convert a Kratos.Matrix to a Numpy array?

Not yet....

[loumalouomega]

This seems to be solved, but I do not understand how I can get the stiffness and mass matrix to Python. Is it documented anywhere?

The LHS is easy, the solvers could be adapted in order to retrieve the K and M in separated way

I just checked and it retrieves in matrix market format (*.mm files with an echo level of 4)

[veiguf]

The problem has surfaced for us again. We would like to access the mass and stiffness matrices of a structural model and get them in Python. How can we do that? Any help would be appreciated!

[RiccardoRossi]

dear @veiguf take a look at

https://github.com/KratosMultiphysics/Kratos/blob/master/applications/StructuralMechanicsApplication/python_scripts/structural_mechanics_custom_scipy_base_solver.py

it provides you essentially what you need in scipy format

and

https://github.com/KratosMultiphysics/Kratos/blob/master/kratos/python_scripts/scipy_conversion_tools.py

to understand how the conversion to scipy happens

[veiguf]

Thanks @RiccardoRossi for pointing out that file, did not know about it... Based on this, we tried a couple things to write a class to extract the system matrices, in vain. TryExtract.zip Should this process be carried out with an arbitrary model or must the project parameters be specifically set for this. Short: Can I build the model with GiD and then use these commands? Is there a more detailed example somewhere where something like this was done?

[RiccardoRossi]

i am not sure i understand you. The idea behind this was to create a mdpa and ProjectParameters as you would for a dynamic analysis, and then run this strategy instead of the normal one.

By looking very quickly in the code of the TryExtract i see that there is a mismatch between "Mass" and "Stiff" which are camel case and "mass" and "stiff" which have no capitals. not sure if this may affect

@manuelmessmer can u comment more on this? you played around with this more than i did

[manuelmessmer]

Hi @veiguf, if I understand it right, you want to set up a model in GiD and then retrieve the system matrices in python, correct? As @RiccardoRossi said, this is exactly what the mentioned structural_mechanics_custom_scipy_base_solver.py does. There is also a test case, which probably answers your questions how it is used. https://github.com/KratosMultiphysics/Kratos/blob/master/applications/StructuralMechanicsApplication/tests/test_custom_scipy_base_solver.py You just need to define the respective solver in the ProjectParameters.json. https://github.com/KratosMultiphysics/Kratos/blob/5566a9cbdc52ed621117d50698a9e88a38f81c1a/applications/StructuralMechanicsApplication/tests/test_custom_scipy_base_solver.py#L33 Then you can also derive the class CustomScipyBaseSolver and override the SolveSolutionStep method and implement there what you want.

[Theshahzeb]

hello, I was assigned this task and to start i ran this test_custom_scipy_base_solver.py test case AS IS but it failed screenshot I am fairly new and I would really appreciate any help. Thanks

[manuelmessmer]

hi, in the screenshot I saw that you called a file called test_compare_scipy_eigen.py. Did you create that file yourself. Seems like your system of equations is singular. If you run python3 test_custom_scipy_solver.py it should work.

[Theshahzeb]

Its the same file actually, replacing props[KratosMultiphysics.COMPUTE_LUMPED_MASS_MATRIX] = False with props[StructuralMechanicsApplication.USE_CONSISTENT_MASS_MATRIX] = True seems to work except these warnings : custom_scipy_test_warnings

[manuelmessmer]

It complains because if a block_builder_and_solver is used the first eigenmodes depend on the applied Dirichlet boundary conditions and will lead to zero eigenvalues. Therefore, in this case it is better to use elimination builder use_block_builder=false under builder_and_solver_settings. The variable in this test is deprecated. It is still block_builder = False. Anyways it should work with the elimination builder (without warning). Btw are you interested in the eigenvalue analysis or do you want do modify the matrices in python for some other reason? Because the eigenvalue analysis is basically only an example to show how this solver works. If you want to do something else, you can override the function SolveSolutionStep in the structural_mechanics_custom_scipy_solver.py.

[manuelmessmer]

I did the changes here: #8355

[Theshahzeb]

Hi @manuelmessmer thank you so much for the clarification and no, I am not interested in the eigenvalues, just trying to get the mass and stiffness matrices from a model to manipulate them in python.

[Theshahzeb]

Hi, I tried to extract mass and stiffness matrices of a model with no boundary conditions using structural_mechanics_custom_scipy_base_solver.py but i get a type error that says SetToZeroMatrix(): incompatible function arguments. Parameters file: ProjectParameters.json. Thanks

[manuelmessmer]

How do you run your analysis? With

simulation = StructuralMechanicsAnalysis(model,parameters)
simulation.Run() 

[Theshahzeb]

I am running structural_mechanics_scipy_base_solver.py (in the GiD model folder) with modified SolveSolutionStep and

with open("ProjectParameters.json", 'r') as parameter_file:
     parameters = KratosMultiphysics.Parameters(parameter_file.read())
model = KratosMultiphysics.Model()
simulation = CustomScipyBaseSolver(model, parameters)
simulation.SolveSolutionStep()

[manuelmessmer]

This will not work, because you need an analysis_stage (In your case StructuralMechanicsAnalysis(model,parameters)) to set up your solvers properly. Then you can start the simulation with Run(). All functions such as InitializeSolutionStep(), SolveSolutionStep() etc.. will then be called in there. In your case the previous steps are not called and the matrices are not constructed yet. Most likely the GetSystemMatrix() was just returning an empty pointer. You don't need to call the CustomScipySolver directly. The analysis_stage will read trough your ProjectParameter.json and invoke the solver listed there.

[Theshahzeb]

hello, i apologize for the delayed response, the analysis_stage solved the problem, reads the files correctly and runs without errors but it does not go into the functions like _MassMatrixComputation & SolveSolutionStep, which means I cannot see any matrix or variable. I really appreciate the help. Thank you

[manuelmessmer]

It might be that your python script is not in the correct folder. Basically you have two options. Either you modify the original base_solver and recompile (it will move it then to the bin folder, from where it is called). Or you keep everything in your new file structural_mechanics_custom_scipy_base_solver2.py, but then you have to specify this file in your ProjectParameter.json. If you do so, it makes sense to derive a new class from the base_solver. Most likely right now you only start the simulation in your new file.

Btw you can put these lines

with open("ProjectParameters.json", 'r') as parameter_file:
    parameters = KratosMultiphysics.Parameters(parameter_file.read())
model = KratosMultiphysics.Model()
# simulation = CustomScipyBaseSolver(model, parameters)
simulation = StructuralMechanicsAnalysis(model, parameters)        
simulation.Run()

into a separate MainKratos.pyfile and start everything from there. Maybe a bit more convenient and cleaner.

[Theshahzeb]

Thank you so much, this makes sense now. I'd go for the 2nd option but i specify my new file structural_mechanics_custom_scipy_base_solver2.py here in the ProjectParameters.json

"solver_type" : "KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_custom_scipy_base_solver",

right?

[manuelmessmer]

Yes, but your file must still be in the correct folder. Otherwise, Kratos will not find it. The other option would be easier, but that's up to you.

[e-dub]

I just implemented a "new" solver, i.e. the "other option" to do this and I can indeed extract the matrices correctly: [new solver].(https://github.com/KratosMultiphysics/Kratos/blob/struct/extract-ssysem-matrices/applications/StructuralMechanicsApplication/python_scripts/structural_mechanics_system_matrix_extraction_solver.py)

I am now trying to think of the best way to get these then into the python environment. Suggestions?

Also I called this thing preliminarily: structural_mechanics_system_matrix_extraction... Suggestion on this?

[manuelmessmer]

I am not sure if I understand correctly. In the SolveSulutionStep() you should have the matrices in python then? Or what is your final goal with the matrices?

[e-dub]

Maybe it is a mistake in my thinking. I would like to have a general extraction method that I can use with any model. What is the easiest way to pass it from the solver back into the "main" Python file, which I use to call the model? ...maybe I am overthinking this though (or underthinking...Friday).

[manuelmessmer]

Well, the solver was "only" designed to solve/ manipulate the matrices with scipy etc.. So you still have the normal Kratos wokrflow and can run nonlinear SolutionLoops. So obviously, you can just pass the matrices from the SolveSolutionStepto any python function. However, to actually return them you would have to override the RunSolutionLoop of the analysis stage. Or call the functions below individually inside your MainKratos (instead of: analysis.Run()). Line 65 calls the SolveSolutionStep of the scipy solver.

https://github.com/KratosMultiphysics/Kratos/blob/c6c48df88b60b0a719d0dc79dab447c7410ccc13/kratos/python_scripts/analysis_stage.py#L43-L68