Failing Example in release

[agiuliodori]

Dear @KlausBSautter @philbucher @loumalouomega , it is Riccardo here writing on behalf of Agustina.

we have a test which is still failing with the beam. It sort of works for the small displacmeent beam, however it fails immediately as soon as we switch to the large displacement one.

https://www.dropbox.com/s/t05ixe69v2jtpq5/Test-beam-metro.gid.zip?dl=0

the test represents a small portion of a stent, with realistic data for the properties of the beams. Data are in mm.

essentially the boundary conditions pull the "tube" in the direction y and fix the displacement of the basis. The interest is essentially in measuring the change of diameter as a function of the relative displacement in y direction.

With the small deformation beam we get a result however the structure becomes bigger in all directions instead of showing a reduction in diameter as we increase the pull. the large dispalcement case fails invariably, even if we use smaller steps. We have been taking some good care in the BCs and we are convinced that they are ok, so in our opinion (Riccardo's in particular :-p) the error should not be in he BCs but rather somewhere in the solution process.

any hint would be very welcome

[agiuliodori]

Better try with this:

https://www.dropbox.com/s/ymamgax8ce5am8i/Test-beam-metro2.gid.zip?dl=0

[KlausBSautter]

With the small deformation beam we get a result however the structure becomes bigger in all directions instead of showing a reduction in diameter as we increase the pull.

This is very strange, please check your settings, I get the following result:

[RiccardoRossi]

@KlausBSautter this is what i would expect...

can u post the settings you employed? (best of all a zip which i can run directly on my computer)

[RiccardoRossi]

@KlausBSautter just to tell that i tried out Agustina's case on my computer in linux and i also miss convergence

[philbucher]

are you on release or master-branch?

also @KlausBSautter probably already used #4478

[KlausBSautter]

@philbucher no I have not used #4478, this only fixed problems with mpcs. @RiccardoRossi I just used the linear beam with a linear analysis. Non-linear does not work for me, neither

[RiccardoRossi]

... so we have a broken nonlinear beam in the release? :-(

the second example of Agustina did work for the linear case as you show...

[KlausBSautter]

roll_up.gid.zip

[KlausBSautter]

snap_through.gid.zip

[KlausBSautter]

2 highly nonlinear cases for which the beam works. Also as mentioned @philbucher , @AndreasWinterstein etc. use it in FSI and other things where it works too. We might have to put in additional work to make it bullet proof

[loumalouomega]

Maybe the problem is on the small beams added for rotation stiffness. @agiuliodori can you switch the short beams with linear beams instead of non-linear beams?

[agiuliodori]

@loumalouomega do you refer to change in the short beams to a beam small displacement element? I tried it and it doesnt work too.

[loumalouomega]

Yes, I meant that.

El mar., 19 mar. 2019 11:32, Agustina notifications@github.com escribió:

@loumalouomega https://github.com/loumalouomega do you refer to change in the short beams to a beam small displacement element? I tried it and it doesnt work too.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/KratosMultiphysics/Kratos/issues/4479#issuecomment-474294850, or mute the thread https://github.com/notifications/unsubscribe-auth/ATEMgNnoSvESej1LuNn_Ib3zoC885MMeks5vYLy7gaJpZM4b6LWv .

[KlausBSautter]

if you do everything small dsip and use a linear analysis it works. I dont know what the problem here is with the large def

[RiccardoRossi]

can u try some tests which involve torsion and that are NOT contained in the XY plane?

[KlausBSautter]

yes wait a second I am doing a case

[RiccardoRossi]

I tried the roll up case rotated by 45 deg. it works for a while but it fails pretty soon

On Tue, Mar 19, 2019 at 2:45 PM Klaus Bernd Sautter < notifications@github.com> wrote:

yes wait a second I am doing a case

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/KratosMultiphysics/Kratos/issues/4479#issuecomment-474376432, or mute the thread https://github.com/notifications/unsubscribe-auth/AHr7EUEFLAbfinhOn49m2kMdQK0CxT4Pks5vYOnjgaJpZM4b6LWv .

--

Riccardo Rossi

PhD, Civil Engineer

member of the Kratos Team: www.cimne.com/kratos

Associate Professor at Universitat Politècnica de Catalunya, BarcelonaTech (UPC)

Full Research Professor at International Center for Numerical Methods in Engineering (CIMNE)

C/ Gran Capità, s/n, Campus Nord UPC, Building C1, First Floor

08034 – Barcelona – Spain – www.cimne.com -

T.(+34) 93 401 56 96 skype: rougered4

http://www.cimne.com/

https://www.facebook.com/cimne http://blog.cimne.com/ http://vimeo.com/cimne http://www.youtube.com/user/CIMNEvideos http://www.linkedin.com/company/cimne https://twitter.com/cimne

Les dades personals contingudes en aquest missatge són tractades amb la finalitat de mantenir el contacte professional entre CIMNE i voste. Podra exercir els drets d'accés, rectificació, cancel·lació i oposició, dirigint-se a cimne@cimne.upc.edu. La utilització de la seva adreça de correu electronic per part de CIMNE queda subjecte a les disposicions de la Llei 34/2002, de Serveis de la Societat de la Informació i el Comerç Electronic.

Imprimiu aquest missatge, només si és estrictament necessari. http://www.cimne.com/

[loumalouomega]

I tried the roll up case rotated by 45 deg. it works for a while but it fails pretty soon

Then the problem is on the rotation matrix?, this is a good new?

[KlausBSautter]

tors.gid.zip @RiccardoRossi

[KlausBSautter]

I tried the roll up case rotated by 45 deg. it works for a while but it fails pretty soon

can you add your test

[KlausBSautter]

I remember that we had this problem already a while ago but I thought we fixed it with the quaternions

[RiccardoRossi]

roll_up_diag.gid.zip

here goes the roll up test at 45deg

[KlausBSautter]

roll_up_diag.gid.zip here goes the roll up test at 45deg

I see, yes that's strange. Especially as the penultimate solution step converges within 3 iterations

[loumalouomega]

Maybe would be interesting to create a c++ test with one beam and one load. We do X rotations of this beam and check that the result is always the same, at least each 45º degrees. That would help to find the error

[agiuliodori]

@RiccardoRossi Here is attached a similar example but with a geometry without the "joints elements". test-without-joints-elements.zip

[KlausBSautter]

@RiccardoRossi I think I found a problem. In the geometric K the internal moments are added whereas I just saw that the book which I used says those should be the external moments at the node.... I find that strange to be true. I removed them and now I can roll up the beam in every possible configuration. That's good news. But for example the "torsional problem" I posted here

tors.gid.zip

now needs more iterations to converge.

@agiuliodori and @RiccardoRossi the examples you are posting here behave very strange. If you change the prescribed displacements to zero. It is still iterating and calculating a residual? Where does this come from? What exactly are you doing here?

With the change in K the element does not crash anymore. But I am not completly happy with these changes as e.g.

tors.gid.zip

needs more iterations

[loumalouomega]

@RiccardoRossi I think I found a problem. In the geometric K the internal moments are added whereas I just saw that the book which I used says those should be the external moments at the node.... I find that strange to be true. I removed them and now I can roll up the beam in every possible configuration. That's good news. But for example the "torsional problem" I posted here

tors.gid.zip

now needs more iterations to converge.

@agiuliodori and @RiccardoRossi the examples you are posting here behave very strange. If you change the prescribed displacements to zero. It is still iterating and calculating a residual? Where does this come from? What exactly are you doing here?

With the change in K the element does not crash anymore. But I am not completly happy with these changes as e.g.

tors.gid.zip

needs more iterations

What do you mean by external?, that are not internal forces?, that must be separated and differentiated? (this remembers me many discussion about separating internal and external forces)

[RiccardoRossi]

Hi @KlausBSautter just to tell that the "end point moments" are defined in equation 5.2 as the "conjugate forces".

To my understanding this means that they are -RHS as obtained for the given displacement and rotations

[KlausBSautter]

yes this is the current state of the element.

[RiccardoRossi]

mmm just try with the opposite sign :-D

[loumalouomega]

mmm just try with the opposite sign :-D

That would be funny

[KlausBSautter]

@RiccardoRossi I'll try tmr. Also the problem posted here https://www.dropbox.com/s/ymamgax8ce5am8i/Test-beam-metro2.gid.zip?dl=0

does iterate for nl-beam but not for the linear beam if no external forces/disp are applied. I have no idea what is happening there and no explanation for that....

[RiccardoRossi]

@KlausBSautter one important comment, the m is the "beam level" reaction to my understanding, NOT the assembled one

[KlausBSautter]

yes looking at a single element right now

[philbucher]

@RiccardoRossi (and @loumalouomega ) Being also involved in this I feel like I have to add my thoughts: in my opinion the problems you are having here are not really the fault of @KlausBSautter

4286 did MAJOR changes to how the convergence-things are checked in the NR-strategy. Before things were working afaik. Esp @AndreasWinterstein, myself and others at the chair have used the beams for large and complex examples without problems.

Now I am not saying that the beam is completely correct implemented, but this came to light only with #4286

Please in the future try to avoid pushing such changes very close to the release, even it seems to be necessary, it created also a lot of problems bcs it had to be done "fast". Same for examples, please try such things also in the "off-season", and not only with having the release next door.

We will speak abt this next week too

[loumalouomega]

@RiccardoRossi (and @loumalouomega ) Being also involved in this I feel like I have to add my thoughts: in my opinion the problems you are having here are not really the fault of @KlausBSautter

4286 did MAJOR changes to how the convergence-things are checked in the NR-strategy. Before things were working afaik. Esp @AndreasWinterstein, myself and others at the chair have used the beams for large and complex examples without problems.

Now I am not saying that the beam is completely correct implemented, but this came to light only with #4286

Please in the future try to avoid pushing such changes very close to the release, even it seems to be necessary, it created also a lot of problems bcs it had to be done "fast". Same for examples, please try such things also in the "off-season", and not only with having the release next door.

We will speak abt this next week too

We don't detect huge bugs just before release on purpose..., I am sorry to be such a trouble. The thing is that the previous implementation could be hiding errors. Don't worry we will solve all this issues on time (I hope)

[RiccardoRossi]

I may agree that it was too big change to be pushed so late in the release. Still i think this was an important bugfix, as it simply prevented convergence in many cases for which convergence was actually met in 1iteration.

Regarding the problem at hand (the nonlinear beam) i believe it is not convergence related. The beam converges nicely until it fails abruptly...and the tests still run fine in the XY plane but fail when simply rotated at 45deg.

We were speaking now in the @KratosMultiphysics/technical-committee to make a bugfix release right after the workshop not to merge latest stuff so late.

[KlausBSautter]

@RiccardoRossi I now went through everything again and really can't find an obvious mistake. BUT I found the problem with the example in this post, why it has a RHS even though I am turning off all forces and prescribed displacement. It seems to be a numeric problem, I tried it with other tools and I can observe the same problem:

I have a vector nx = (0.3678038875014457698142678055,-0.8716178212545082626405701376,-0.3240411609816999938260551062 and a vector n = (0.3678038875014457698142678055,-0.8716178212545082626405701376,-0.3240411609816999938260551062)

If there is no displacement n=nx as it is calculated correctly by the beam.

if I do the following calculation (forget about ny,nz):

You expect the new nx to be equal to the old nx, which is not happening, the result is nx = (0.3678038875014458253254190367,-0.8716178212545084846851750626,-0.3240411609817001048483575687)

which is different after some digits. This difference comes into play when calculating the new internal forces. Especially for small example as the one in this post this will lead to internal forces >1e-3 and will make this system not converge. I have a feeling that this numerical innacurray might also lead to problems when the internal moments are calculated for a beam structure which is rolled up in a different setup than [1,0,0] || [0,1,0] || [0,0,1].

The corresponding lines in the beam are: https://github.com/KratosMultiphysics/Kratos/blob/f65f4082d79c1b4cdd6f4feb35c95c83c8288d84/applications/StructuralMechanicsApplication/custom_elements/cr_beam_element_3D2N.cpp#L856-L874

I tested this with numpy and a high precision and do observe the same behaviour. Do you think using sth. else but outer_prod and prod solves this?

[RiccardoRossi]

well, the implementation is not optimal, but it is not what we are speaking about a good implementation (much faster btw) would be

double tmp = inner_prod(Bisectrix,n_xyz);
noalias(n_xyz) -= 2.0*tmp*Bisectrix;

in any case why do you expect it to be equal? It might be so in the xy plane, but i don't see it obvious in 3D

[KlausBSautter]

It must be equal because no displacement is present

[KlausBSautter]

so this is a big poblem because exactly this causes really small examples (as this one) to produce inner forces without external influence

[RiccardoRossi]

i mean, the vectors you are putting here have norm 1, so i don't think it is a numerical error

import numpy
nx = numpy.array([0.3678038875014457698142678055,-0.8716178212545082626405701376,-0.3240411609816999938260551062])
n = numpy.array([0.3678038875014457698142678055,-0.8716178212545082626405701376,-0.3240411609816999938260551062])
y = numpy.array([-nx[0],nx[1],nx[2]])
print(n-nx)

tmp = numpy.dot(n,y)
print("dot(n,y)=",tmp)
a = y - 2.0*tmp*n 
print("res = ",a)

outputs

[0. 0. 0.]
dot(n,y)= 0.7294406006776478
res =  [-0.90438606  0.39996903  0.1486964 ]

which is VERY different from nx...

[KlausBSautter]

import numpy as np
np.set_printoptions(precision=50)
n = np.array([0.3678038875014457698142678055,-0.8716178212545082626405701376,-0.3240411609816999938260551062])
I = np.eye(3)
n_outer = np.outer(n,n)
a1 = I-2.0*n_outer
a2 = a1@(-n)

n = [ 0.36780388750144577 -0.8716178212545083  -0.3240411609817    ]
a2 = [ 0.3678038875014458 -0.8716178212545085 -0.3240411609817001]
cross(a2,n) = [0.0000000000000000e+00 1.3877787807814457e-17 0.0000000000000000e+00]

[KlausBSautter]

I don't know exatly what you are doing, here you see that a1 should be equal to n but is not excatly the same

[RiccardoRossi]

oh, it is multiplied by -n (you wrote (-nx,ny,nz) in the previosu post

[KlausBSautter]

yes it is -nx, I just use n=nx here

[RiccardoRossi]

import numpy
nx = numpy.array([0.3678038875014457698142678055,-0.8716178212545082626405701376,-0.3240411609816999938260551062])
n = numpy.array([0.3678038875014457698142678055,-0.8716178212545082626405701376,-0.3240411609816999938260551062])
y = -nx
print(n-nx)

tmp = numpy.dot(n,y)
print("dot(n,y)=",tmp)
a = y - 2.0*tmp*n 
print("res = ",a)
print("res-nx",a-nx)

gives

[0. 0. 0.] dot(n,y)= -1.0000000000000002 res = [ 0.36780389 -0.87161782 -0.32404116] res-nx [ 1.11022302e-16 -4.44089210e-16 -1.11022302e-16]

[KlausBSautter]

import numpy as np
np.set_printoptions(precision=50)
n = np.array([0.3678038875014457698142678055,-0.8716178212545082626405701376,-0.3240411609816999938260551062])
nx = np.array([0.3678038875014457698142678055,-0.8716178212545082626405701376,-0.3240411609816999938260551062])
I = np.eye(3)
n_outer = np.outer(n,n)
a1 = I-2.0*n_outer
a2 = a1@(-nx)
print(n)
print(nx)
print(a2)
print(np.cross(n,a2))

n=[ 0.36780388750144577 -0.8716178212545083  -0.3240411609817    ]
nx=[ 0.36780388750144577 -0.8716178212545083  -0.3240411609817    ]
a2=[ 0.3678038875014458 -0.8716178212545085 -0.3240411609817001]
cross(n,a2) =[0.0000000000000000e+00 1.3877787807814457e-17 0.0000000000000000e+00]

[KlausBSautter]

res-nx [ 1.11022302e-16 -4.44089210e-16 -1.11022302e-16]

yes should be 0

[KlausBSautter]

increment_deformation : [12](0,0,0,0,0,0,0,0,0,0,0,0)
<<
delta_x : [3](0.367804,-0.871618,-0.324041)
0.3678038875014457698142678055 | -0.8716178212545082626405701376 | -0.3240411609816999938260551062
rotated_nx0 : [3](0.3678038875014457698142678055,-0.8716178212545082626405701376,-0.3240411609816999938260551062)
0.3678038875014457698142678055 | -0.8716178212545082626405701376 | -0.3240411609816999938260551062
<<
rotated_nx0 : [3](0.3678038875014458253254190367,-0.8716178212545084846851750626,-0.3240411609817001048483575687)
Bisectrix : [3](0.3678038875014457698142678055,-0.8716178212545082626405701376,-0.3240411609816999938260551062)
temp_vector : [3](0,-1.387778780781445675529539585e-17,0)
deformation_modes_total_v : [6](0,0,0,0,-2.158170144853027542259783778e-17,1.657279071734585126883546443e-17)
element_forces_t : [6](0,0,0,0,-3.374323339415775198282115843e-05,2.591174502629473282784740296e-05)

you can see that the temp_vector is the result of a a crossproduct between bisectrix and rotated_nxo which is not the same anymore due to the numerical errror. This temp_vector is used to calcualte internal forces element_forces_t with the help of the element stiffnesses