Kratos: proposal of a new option for the EigensolverStrategy

[RiccardoRossi]

@armingeiser @oberbichler as commented in #6099 i am trying to provide a way to assemble K, M and eventually D and make them available as scipy matrices.

as far as i understand the computation of K and M, with the related dofs is done in the EigenSolverStrategy.

I was wondering if you would agree about an additional option for the strategy on the line of
"only_matrix_assemble"

this would be false by default, but if one sets it to true than the eigenvalue simulation and the writing of the eigenvalues/vectors to kratos would be skipped. However one could retrieve the K and M matrices and the dofs, convert them to scipy and work with them.

what do you think?

[armingeiser]

For me that would be fine.

However if at some point you need D as well, i don't think the eigensolver_strategy is the right place for that.
I am not that much familiar with the strategies, but would it be possible to create a simple function for the assembly of K, M, (D) respectively, that can be used from within the eigen_solver strategy (and potentially also other strategies). Assembly of the matrices should always be the same right?
Or could you get the matrices from the BuilderAndSolver directly?

That would be cleaner/clearer IMO for a general way to get those matrices for a system in python.
If the strategies modify them somehow after they get them from the builder and solver it would not work.

[armingeiser]

About the scipy interface:

Are you considering the Eigen/pybind11/numpy/scipy connection? As scipy is already a dependency it would probably not be that bad to add the Eigen dependency as well (for this feature) and make use of the already existing interface through pybind

[philbucher]

+++1 for @armingeiser
One more reason to put Eigen in the core

[RiccardoRossi]

i am not against having eigen in the core, and @philbucher opened an issue about that so i would not discuss it here but rather in that issue.

My remark here is that scipy IS NOT a dependency (also numpy is not, as of now).
The outstanding blocker (see #2165 ) is that it is highly problematic to distribute them if we have our own python wrapper (which we need to distribute kratos together with GiD). The recent change of installation mechanism made the kratos installation much more pythonic, which in turn should use nuitka or similar to package a launcher. nevertheless i don't think doing it in a robust way is trivial.

Regarding the specific issue, i can see your point. The thing is that in order to correctly compute K and M (as well as D) you need to do all the Initialization steps that are done in the EigenStrategy.
For example elements and conditions need to be initialized, processes need to be applied, etc etc.

At that point the function you are speaking about is effectively rewriting everything in the strategy aside of the SolveSolutionStep.
One option would be to derive a CustomEigenSolvingStrategy from the EigenSolverStrategy and redefine the function SolveSolutionStep from the python library, to do something like

~~~py
def _MassMatrixComputation(self):
... here we will need to resize and initialize M,K, D ...

      self.model_part.ProcessInfo[BUILD_LEVEL] = 1 #this is pretty ugly btw   
      self.eigen_strategy.GetBuilderAndSolver().BuildLHS(self.eigen_scheme, self.model_part, M)
      mass_matrix = KratosMultiphysics.scipy_conversion_tools.to_csr(M)
      M.Clear()
     return mass_matrix


def _StiffnessMatrixComputation(self): 
      ... ResizeAndInitialize(....K ... )
      self.model_part.ProcessInfo[BUILD_LEVEL] = 2 #this is pretty ugly btw   
      self.eigen_strategy.GetBuilderAndSolver().BuildLHS(self.eigen_scheme, self.model_part, K)
      stiffness_matrix = KratosMultiphysics.scipy_conversion_tools.to_csr(K)
      K.Clear()
     return stiffness_matrix

def _DampingMatrix(self): 
     ...same...

def SolveSolutionStep(self):
     M = self._MassMatrix()
     K = self._StiffnessMatrix() 

     dofs = self.strategy.GetDofsArray()

     #here operate as you like with M and K (which are here scipy matrices)

~~~

anyway, i will think about how to obtain this clean

[RiccardoRossi]

well, a prototype would look much similar to the following:

~~~py
from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7

Importing the Kratos Library

import KratosMultiphysics

Import applications

import KratosMultiphysics.StructuralMechanicsApplication as StructuralMechanicsApplication

Import base class file

from KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_solver import MechanicalSolver

from KratosMultiphysics import eigen_solver_factory

def CreateSolver(main_model_part, custom_settings):
return CustomScipySolver(main_model_part, custom_settings)

class CustomScipySolver(MechanicalSolver):
"""The structural mechanics eigen solver.

This class creates the mechanical solvers for eigenvalue analysis.

See structural_mechanics_solver.py for more information.
"""
def __init__(self, main_model_part, custom_settings):
    # Construct the base solver.
    super(CustomScipySolver, self).__init__(main_model_part, custom_settings)
    KratosMultiphysics.Logger.PrintInfo("::[CustomScipySolver]:: ", "Construction finished")

@classmethod
def GetDefaultSettings(cls):
    this_defaults = KratosMultiphysics.Parameters("""{
        "scheme_type"         : "dynamic"
    }""")
    this_defaults.AddMissingParameters(super(CustomScipySolver, cls).GetDefaultSettings())
    return this_defaults

#### Private functions ####

def _create_solution_scheme(self):
    """Create the scheme for the eigenvalue problem.

    The scheme determines the left- and right-hand side matrices in the
    generalized eigenvalue problem.
    """
    scheme_type = self.settings["scheme_type"].GetString()
    if scheme_type == "dynamic":
        solution_scheme = StructuralMechanicsApplication.CustomScipySolverDynamicScheme()
    else: # here e.g. a stability scheme could be added
        err_msg =  "The requested scheme type \"" + scheme_type + "\" is not available!\n"
        err_msg += "Available options are: \"dynamic\""
        raise Exception(err_msg)

    return solution_scheme

def _create_linear_solver(self):
    return null

def _create_mechanical_solution_strategy(self):
    eigen_scheme = self.get_solution_scheme() # The scheme defines the matrices of the eigenvalue problem.
    builder_and_solver = self.get_builder_and_solver() # The CustomScipySolver is created here.
    computing_model_part = self.GetComputingModelPart()

    return KratosMultiphysics.LinearStrategy(computing_model_part,
                                            eigen_scheme,
                                            builder_and_solver,
                                            self.settings["compute_modal_decomposition"].GetBool())

def _MassMatrixComputation(self): 
    self.GetComputingModelPart().ProcessInfo.SetValue(KratosMultiphysics.BUILD_LEVEL,1)
    scheme = self.get_mechanical_solution_strategy().GetScheme() 

    b = self.get_mechanical_solution_strategy().GetSystemVector()
    self.space.SetToZero(b)

    aux = self.get_mechanical_solution_strategy().GetSystemMatrix() 
    self.space.SetToZero(aux)

    builder_and_solver = self.get_mechanical_solution_strategy()
    builder_and_solver.Build(scheme, self.GetComputingModelPart(), aux, b)
    builder_and_solver.ApplyConstraints(scheme, self.GetComputingModelPart(), aux, b)

    #guess we need to define an utility function to applies nicely dirichlet conditions
    builder_and_solver.ApplyDirichletConditions(...)

    M = KratosMultiphysics.scipy_conversion_tools.to_csr(aux)
    return M

def _StiffnessMatrixComputation(self): 
    self.GetComputingModelPart().ProcessInfo.SetValue(KratosMultiphysics.BUILD_LEVEL,2)
    scheme = self.get_mechanical_solution_strategy().GetScheme() 

    b = self.get_mechanical_solution_strategy().GetSystemVector()
    self.space.SetToZero(b)

    aux = self.get_mechanical_solution_strategy().GetSystemMatrix() 
    self.space.SetToZero(aux)

    builder_and_solver = self.get_mechanical_solution_strategy()
    builder_and_solver.Build(scheme, self.GetComputingModelPart(), aux, b)
    builder_and_solver.ApplyConstraints(scheme, self.GetComputingModelPart(), aux, b)

    #guess we need to define an utility function to applies nicely dirichlet conditions
    builder_and_solver.ApplyDirichletConditions(...)

    K = KratosMultiphysics.scipy_conversion_tools.to_csr(aux)
    return K

def SolveSolutionStep(self):
    #obtain scipy matrices
    M = self._MassMatrixComputation()
    K = self._StiffnessMatrixComputation()

    dofs = self.get_mechanical_solution_strategy().GetDofSet()

    #...here one implement stuff with scipy

    return true #converged

~~~

[armingeiser]

At that point the function you are speaking about is effectively rewriting everything in the strategy aside of the SolveSolutionStep.

I see your point. And this would need to be done for every strategy one potentially wants to have the matrices of...
Then maybe a flag in the strategy that skips the solution is not too bad. And it could be easily integrated in other strategies as well.

Your proposal about the CustomScipySolver does more then just providing the matrices. If anyway what you want to do is solution of eigenvalue problems with scipy, its good, if its only getting the matrices i think it is overkill.

[qaumann]

@RiccardoRossi on the mor-application brach, there is system_matrix_builder_and_solver.hpp which assembles the system matrices (K,C,M). @adityaghantasala implemented it and I (will) use it for my model order reduction methods, perhaps this is similar to what you need?

[RiccardoRossi]

@qaumann what is the difference between that builder and solver and the one
in the core?

On Thu, Dec 19, 2019, 1:47 PM qaumann [e-mail] wrote:

@RiccardoRossi https://github.com/RiccardoRossi on the mor-application
brach, there is system_matrix_builder_and_solver.hpp which assembles the
system matrices (K,C,M). @adityaghantasala
https://github.com/adityaghantasala implemented it and I (will) use it
for my model order reduction methods, perhaps this is similar to what you
need?

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/KratosMultiphysics/Kratos/issues/6111?email_source=notifications&email_token=AB5PWEOGRN7SSROTVJT7N5DQZNUO5A5CNFSM4J35XHG2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHJP56Q#issuecomment-567475962,
or unsubscribe
https://github.com/notifications/unsubscribe-auth/AB5PWEJADBHW2WZOVSHBLWTQZNUO5ANCNFSM4J35XHGQ
.

[RiccardoRossi]

i ll start by saying that i was originally very much in favor (aside of
licensing doubts) of switching to eigen, I also agree that ublas is a dead
end.
also, i have been reading about both eigen and blaze.

long sotry shirt, I am currently (much) less clear that it is a good idea
than i was some time ago.

here go a few detais that i see as problematic (for a deep integration,
more on this later):
1 - eigen is essentially a c++03 lib. examples of this is the lack of usage
of initializer lists (according to docs), problematic use of "auto" and, as
i understand, lack of use of modern alignement mechanisms (c++11-17).
i wonder if we woukd not be jumping onto an already old lib
2 - i see as a blocker that eigen disallows pass by value. independently on
considering if it is a good idea or not (personally, i do not see anything
wrong about doing it, if what you want is a copy), it is a feature
currently used through the code. without compiler support i do not see any
way of fixing this systematically.
3 - i am not clear about openmp in eigen. we would need openmp for sparse
matrices but not for dense ones (which for us are pretty tiny). as i
understand one can either activate it for all eigen or completely disallow
it.
4 - we woukd need a compatibility layer the same way we did for Amatrix. i
am not sure of the implications/feasibility of it in the case of Eigen.
5 - i think that for maximal performance (and strictly optionally) eigen
uses internally blas for some kernels. i hate the idea of having the core
to depend on blas/lapack.

having said this, and without entering in the pybind/numpy/discussion for
now, i think we have a few open possibilities with different level of
intrusiveness:

1 - (minimally intrusive) leave it as now but packaging eigen by default.
(would need to remain optional, and everything should work without, as now)
2 - use eigen only for sparse linear algebra. this would imply changing
builder and solver, etc, but not all the rest of kratos. on the bright side
eigen sparse would be always available. Question: is there any alternative
aside of eigen?
If 2 is the decision in my opinion we should explicitly disallow using
eigen for dense matrices to avoid it being mixed with the library we use
for tiny matrices. We would then need our own (or another lib) for tiny
matrices.
3 - using eigen for all (or an other lib, say blaze). if that is the
decision, than we should also design a migration strategy, ideally with
compiler support. (runtime segfaults for correct code are not
admissible...). also i would like my comments above to be answered...

finally, and about numpy/scipy if anyone has an idea of how to write a
"runkratos" that supports it on a clean computer, proposals are very
welcome.

for now my only idea is to have a "runkratos.py" which looks something like

from numpy import *
from scipy import *

exec( ... mainkratos.py... )

and to hope that nuitka (or cxfreeze) will correctly package numpy and
scipy dependencies, including blas and lapack (...shiver...)

On Thu, Dec 19, 2019, 1:47 PM qaumann [e-mail] wrote:

@RiccardoRossi https://github.com/RiccardoRossi on the mor-application
brach, there is system_matrix_builder_and_solver.hpp which assembles the
system matrices (K,C,M). @adityaghantasala
https://github.com/adityaghantasala implemented it and I (will) use it
for my model order reduction methods, perhaps this is similar to what you
need?

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/KratosMultiphysics/Kratos/issues/6111?email_source=notifications&email_token=AB5PWEOGRN7SSROTVJT7N5DQZNUO5A5CNFSM4J35XHG2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHJP56Q#issuecomment-567475962,
or unsubscribe
https://github.com/notifications/unsubscribe-auth/AB5PWEJADBHW2WZOVSHBLWTQZNUO5ANCNFSM4J35XHGQ
.

[RiccardoRossi]

oops, i oosted this commebt in the wrong issue. i will also post it in the
correct one

On Sat, Dec 21, 2019, 3:02 PM Riccardo Rossi [e-mail] wrote:

i ll start by saying that i was originally very much in favor (aside of
licensing doubts) of switching to eigen, I also agree that ublas is a dead
end.
also, i have been reading about both eigen and blaze.

long sotry shirt, I am currently (much) less clear that it is a good idea
than i was some time ago.

here go a few detais that i see as problematic (for a deep integration,
more on this later):
1 - eigen is essentially a c++03 lib. examples of this is the lack of
usage of initializer lists (according to docs), problematic use of "auto"
and, as i understand, lack of use of modern alignement mechanisms
(c++11-17).
i wonder if we woukd not be jumping onto an already old lib
2 - i see as a blocker that eigen disallows pass by value. independently
on considering if it is a good idea or not (personally, i do not see
anything wrong about doing it, if what you want is a copy), it is a feature
currently used through the code. without compiler support i do not see any
way of fixing this systematically.
3 - i am not clear about openmp in eigen. we would need openmp for sparse
matrices but not for dense ones (which for us are pretty tiny). as i
understand one can either activate it for all eigen or completely disallow
it.
4 - we woukd need a compatibility layer the same way we did for Amatrix. i
am not sure of the implications/feasibility of it in the case of Eigen.
5 - i think that for maximal performance (and strictly optionally) eigen
uses internally blas for some kernels. i hate the idea of having the core
to depend on blas/lapack.

having said this, and without entering in the pybind/numpy/discussion for
now, i think we have a few open possibilities with different level of
intrusiveness:

1 - (minimally intrusive) leave it as now but packaging eigen by default.
(would need to remain optional, and everything should work without, as now)
2 - use eigen only for sparse linear algebra. this would imply changing
builder and solver, etc, but not all the rest of kratos. on the bright side
eigen sparse would be always available. Question: is there any alternative
aside of eigen?
If 2 is the decision in my opinion we should explicitly disallow using
eigen for dense matrices to avoid it being mixed with the library we use
for tiny matrices. We would then need our own (or another lib) for tiny
matrices.
3 - using eigen for all (or an other lib, say blaze). if that is the
decision, than we should also design a migration strategy, ideally with
compiler support. (runtime segfaults for correct code are not
admissible...). also i would like my comments above to be answered...

finally, and about numpy/scipy if anyone has an idea of how to write a
"runkratos" that supports it on a clean computer, proposals are very
welcome.

for now my only idea is to have a "runkratos.py" which looks something like

from numpy import *
from scipy import *

exec( ... mainkratos.py... )

and to hope that nuitka (or cxfreeze) will correctly package numpy and
scipy dependencies, including blas and lapack (...shiver...)

On Thu, Dec 19, 2019, 1:47 PM qaumann [e-mail] wrote:

@RiccardoRossi https://github.com/RiccardoRossi on the mor-application
brach, there is system_matrix_builder_and_solver.hpp which assembles the
system matrices (K,C,M). @adityaghantasala
https://github.com/adityaghantasala implemented it and I (will) use it
for my model order reduction methods, perhaps this is similar to what you
need?

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/KratosMultiphysics/Kratos/issues/6111?email_source=notifications&email_token=AB5PWEOGRN7SSROTVJT7N5DQZNUO5A5CNFSM4J35XHG2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHJP56Q#issuecomment-567475962,
or unsubscribe
https://github.com/notifications/unsubscribe-auth/AB5PWEJADBHW2WZOVSHBLWTQZNUO5ANCNFSM4J35XHGQ
.

[qaumann]

@qaumann what is the difference between that builder and solver and the one in the core?

It is basically just a builder without solver that assembles the stiffness, mass, and damping matrices. Mass and damping are not assembled in the core residualbased_block_builder_and_solver.

[philbucher]

I think we can close this after #6250
@RiccardoRossi ok for you?