Fractional Matrix Power

NGSolve
1

Hello,

I am currently concerned with some hyperelastic problems, where i want to form some invariants containing arbitrary fractional powers of the right Cauchy Green deformation tensor.

Unfortunately i have found no built in tools which are able to do this.
Therefore i have written some code that allows for the spectral representation of some arbitrary symmetric second order tensor.

However, if i want to do some postprocessing or do my own linearization as well using Assemble it takes forever to assemble the bilinearform.
In attachment you will find the code.

Spectral_Representation.py (5.6 KB)

How can i circumvent this problem or is there an easier and more efficient way to compute and utilize fractional powers of some second order tensor in the Lagrangian?

Thanks in advance and best regards

Philipp

2

Hi Philipp,

there is the possibility to compute the eigenvalues and eigenvectors of a CoefficientFunction using

eigen_vecs_vals = mat.Eig()

It returns a (d+1)d vector, where in the first d*d entries the eigenvectors are stored and in the last d entries the eigenvalues.

This should be faster than computing everything in Python.

Best,
Michael

3

I agree with Michael to go with the built-in Eig - function.

For curiosity - did you look at the expression tree, and did you try with Compile ?
like here:
https://docu.ngsolve.org/latest/i-tutorials/unit-1.2-coefficient/coefficientfunction.html

Joachim

4

Hello Michael & Joachim,

thanks for your fast reply. I am aware of the built-in Eig - function.
However when i define the Lagrangian using a simple Neo Hooke model, where i use the first principle invariant of the right Cauchy Green tensor formulated in terms of its eigenvalues using the built-in Eig - function there seems to be a problem with automatic differentiation and non-distinct eigenvalues.

Since solving a simple problem already gives NaN in the first increment.

@Joachim: No, i did not try with Compile. It seems a little weird to me since AssembleLinearization works perfectly fine with my python code and doing my own Linearization with Assemble or Interpolate 1.PK stress Tensor to some MatrixL2 space does not work.

Best,
Philipp

5

Ok, then let us try where we can get with your Python code.

Can you post a MWE with AssembleLinearization, and a Newton solver ?

Should your own linearization be the exact linearization, what NGSolve also computes ? Or is it some inexact, cheaper approximation ?

What do you mean with interpolating the 1. PK ? Calling Interpolate, maybe to some HDiv**3 space ? Can you give a reference for that method ?

Please just try with .Compile of your expression. Default is ‘soft’ compile, it’s just internal linearization and common sub-expression elimination of the complicated tree. Only calling .Compile(realcompile=True) generates C++ code and calls the compiler.

Joachim

6

Hello Joachim,

thanks again for your fast reply.
I will be on easter holiday and out of office until end of next week.

When i am back i will prepare a MWE in a jupyter-notebook, which illustrates the problem(s).

For postprocessing i called interpolate to some MatrixL2 space.

Ultimate goal is to use ngsolve for some large scale 3D multiphysics (rate-independent dissipative magneto-mechanics etc) [50mio+ DoFs] simulations, where i want to utilize some spectral-type invariants in my consitutive model.

Best and happy easter holidays

Philipp

7

Hello Joachim & Michael,

in attachment you find a MWE that should illustrate the aforementioned problems:
MWE.ipynb (5.6 KB)

  1. There seems to be a problem when formulating the Lagrangian using Eigenvalues of a CF via built-in Eig-function.

  2. Assembling the bilinear form of the given hyperelastic problem using my python code for spectral representation & Assemble() takes forever. The same holds for calling Interpolate for e.g. the 1.PK stress to some MatrixL2 space.

Best

Philipp

8

Hi Philipp,

the NewtonMinimization is below a second, right ?
The interpolation of _PK1 takes forever.
You can significantly improve it by (then it takes 9 sec for me):

PK1_pp.Interpolate(_PK1.Compile())

To get an idea what’s going on you can do:

print (NeoHooke_FP(_F))

which gives about 90k lines, printing all nodes of the expression tree recursively.

You get information of the compiled expression with:

ngsglobals.msg_level = 10
NeoHooke_FP(_F).Compile()

which gives about 150 evaluation steps.

The expression tree of the derivative has 6M nodes, after compilation about 20k steps.

Either use command-line Python, or pipe the stdout to a file, otherwise jupyter-notebook gets stuck.

Maybe some optimizations of the differentiation are missing. We are happy to help you debugging the trees

Joachim

9

Hi Joachim,

thanks for the heads-up. I will have a look at it.
However, for a more efficient implementation it seems wiser to rely on the built-in Eig()-function.

But there seems to be an inherent problem with the Autodiff-functionality in combination with built-in Eig()-function (as also pointed out in my MWE), since it gives wrong results.

Do you have an idea what causes this?

I can also define:

def Matrix_Power(C, n):
    eig_C = C.Eig()
    lbda_1 =  eig_C[9]
    lbda_2 =  eig_C[10]
    lbda_3 =  eig_C[11]
    N_1 = CF( (eig_C[0], eig_C[1],eig_C[2] ) )
    N_2 = CF( (eig_C[3], eig_C[4],eig_C[5] ) )
    N_3 = CF( (eig_C[6], eig_C[7],eig_C[8] ) )
    return lbda_1**n * OuterProduct(N_1, N_1) +  lbda_2**n * OuterProduct(N_2, N_2) + lbda_3**n * OuterProduct(N_3, N_3)

However, when i differentiate expressions involving the built-in Eig()-function the results are always incorrect.

Best

Philipp

10

it is always good to have more options …

with

_PK1 = (NeoHooke_FP(_F)).Diff(_F).Compile()
PK1_pp.Interpolate(_PK1)

I am now down to 0.5 sec. Can you confirm we get the result you expect ?

There was some optimiztion left out, which is now fixed (Diff of atan2 used partial derivative fallback). Before it took 10 sec.

Joachim

11

you are absolutely right and yes, the code runs faster now.

However, there is still the issue with the Autodiff of functions involving the built-in Eig()-function.
(Second option “Flag=0” in my MWE or alternatively the Matrix_Power function that I defined before)

It unfortunately yields false expressions, when I apply automatic differentiation to them.

Best

Philipp

12

Hello again,

just to illustrate the major problem – that autodiff of a CF involving built-in Eig()-function yields completely false results – i attached a minimum working example.
MWE_Differentiation.ipynb (6.5 KB)

To this end i computed some error norms that involves the first principle invariant of the right Cauchy Green Tensor computed in different ways:

## Some error norms ##


_F = F(gfsol).MakeVariable()

_C= _F.trans * _F
_eig_C = _C.Eig()


_I1= Trace(_C)
_I1_MP= Trace(Matrix_Power(_C, 1.0))                 # equivalent to _eig_C[9] + _eig_C[10] + _eig_C[11]
_I1_TP= Trace(Tensor_Power(_C, 1.0)) 


L2_I1_MP= Integrate(sqrt((_I1-_I1_MP)**2), mesh)
L2_I1_TP= Integrate(sqrt((_I1-_I1_TP)**2), mesh)


print(L2_I1_MP)
print(L2_I1_TP)

# Matrix_Power and Tensor_Power yield same results #


l2_I1_F_MP= sqrt(InnerProduct((_I1.Diff(_F)-_I1_MP.Diff(_F)), (_I1.Diff(_F)-_I1_MP.Diff(_F)))).Compile()
l2_I1_F_TP= sqrt(InnerProduct((_I1.Diff(_F)-_I1_TP.Diff(_F)), (_I1.Diff(_F)-_I1_TP.Diff(_F)))).Compile()

L2_I1_F_MP= Integrate(l2_I1_F_MP, mesh)
L2_I1_F_TP= Integrate(l2_I1_F_TP, mesh)

print(L2_I1_F_MP)
print(L2_I1_F_TP)

# Huge error when applying autodiff to built-in Eig() function #

This shows a huge error when i differentiate a CF/GF involving built-in Eig()-function

Hope you can give some help since this should clearly be not the case.

Best

Philipp

13

differntiation auf Eig is not implemented. CF.Eig silently returned wrong results - in the latest version on github it throws an not-available error.

If you want to implement these functions we can give some support.

A starting point is here:

Joachim

14

Thanks Joachim for the clarification. I will have a look at it, when I find the time to do it.

However, differentiation of my own spectral representation code works fine. It is just not the
most performant implementation.

Philipp