Periodic - servo-control using stiffness - simulation loop

1/ This equation is formally correct, I don't see how a gain could be
introduced. Actually, I'd keep calling 1/mass the "inverse inertia"
rather than "gain", since it explains the two methods: one is static,
the other is dynamic. It is better to have two methods, since different
users will prefer one or the other. The only problem in dry granular
materials is that there is no simple derivation of stiffness (unlike
cemented materials with 12 interactions per body or more => no rolling).

2/ The comment is about this cycle :
resetForces->integrateAndUpdate->contactLaw->PeriCompressor->NewtonsLaw.
The simplest way to understand what id does is to rewrite it differently
(it's a cycle, so I can start anywhere) :
- compute forces at time "t" (contactLaw)
- define a new velGrad on the basis of computed forces
- integrate bodies motion based on forces at time "t" and corresponding
velGrad
- integrate cell deformation using the same velGrad (we reseted forces
in between, it doesn't matter)

It should be correct if I don't miss something: forces(t) are used to
define velGrad(t+dt/2) and body.vel(t+dt/2). I was not 100% sure when I
wrote the code, so I added this note. The "else" would be using forces
at time t to define velGrad at t+3dt/2.

Bruno

Hi Bruno,

thanks for answer.

On 30 May 2011 19:01, Chareyre <email address hidden> wrote:

> Your question #159394 on Yade changed:
> https://answers.launchpad.net/yade/+question/159394
>
> Status: Open => Answered
>
> Chareyre proposed the following answer:
> 1/ This equation is formally correct, I don't see how a gain could be
> introduced.

This equation is very similar to the servo-control we use with rigid
boundaries; with the difference that there we also have a gain parameter
(let me call it that way, for me it does not matter if it is split into two
parameters other than one). So, why here the formula is correct but with
rigid walls we need the gain parameter?

> Actually, I'd keep calling 1/mass the "inverse inertia"
> rather than "gain", since it explains the two methods: one is static,
> the other is dynamic. It is better to have two methods, since different
> users will prefer one or the other.

It is not good to have two methods if one of them does not work (at least to
my experience, but I would be happy to hear from others that I am wrong). Am
referring here to what you call "static" method (in PeriIso file). If you
ever try to use it, you will see that it does not properly work. This is why
my question about gain parameter.
Mass is not mass, given its current unit of measure which can be derived
from the equation. As I see it, it is really a gain parameter (and I still
do not see the advantage to use 1/mass when this could simply be called
gain, unless you have a method which estimates the mass parameter from a
physical point of view, do you?).

> The only problem in dry granular
> materials is that there is no simple derivation of stiffness (unlike
> cemented materials with 12 interactions per body or more => no rolling).
>
Ok.

>
> 2/ The comment is about this cycle :
> resetForces->integrateAndUpdate->contactLaw->PeriCompressor->NewtonsLaw.
> The simplest way to understand what id does is to rewrite it differently
> (it's a cycle, so I can start anywhere) :
> - compute forces at time "t" (contactLaw)
> - define a new velGrad on the basis of computed forces
> - integrate bodies motion based on forces at time "t" and corresponding
> velGrad
> - integrate cell deformation using the same velGrad (we reseted forces
> in between, it doesn't matter)
>
> It should be correct if I don't miss something: forces(t) are used to
> define velGrad(t+dt/2) and body.vel(t+dt/2). I was not 100% sure when I
> wrote the code, so I added this note. The "else" would be using forces
> at time t to define velGrad at t+3dt/2.
>
Ok, thanks a lot.

Chiara

>
> Bruno
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/159394/+confirm?answer_id=0
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/159394
>
> You received this question notification because you asked the question.
>

> This equation is very similar to the servo-control we use with rigid
> boundaries; with the difference that there we also have a gain parameter
> (let me call it that way, for me it does not matter if it is split into two
> parameters other than one). So, why here the formula is correct but with
> rigid walls we need the gain parameter?
For rigid walls (or periodic && !dynCell), it reads
displacement=stress_error/stiffness in both cases (a bit more
sophisticated in fact, but this is the basic idea).
For dynCell, it reads acceleration=stress_error/mass (idem).
I don't see what you mean by "we need a gain". Since a gain is any
factor that enter a servo-control, there is basically a gain in all
three cases.

>
> It is not good to have two methods if one of them does not work (at least to
> my experience, but I would be happy to hear from others that I am wrong). Am
> referring here to what you call "static" method (in PeriIso file). If you
> ever try to use it, you will see that it does not properly work. This is why
> my question about gain parameter.

I think Vaclav and Jan used it successfully, with a different
constitutive models than yours.
I know Jan also tried to adjust the stiffness during iterations, instead
of just predicting. So in principle it could work for many constitutive
laws.
I can't comment more personally, I never tried.
> Mass is not mass, given its current unit of measure which can be derived
> from the equation. As I see it, it is really a gain parameter (and I still
> do not see the advantage to use 1/mass when this could simply be called
> gain, unless you have a method which estimates the mass parameter from a
> physical point of view, do you?).
No, I don't have a clear derivation yet. I'm only guessing it needs to
include lengths somewhere. It's an opened question for me (and if I read
this book correctly, it's not explained either
http://www.iste.co.uk/index.php?f=a&ACTION=View&id=384). I suggest you
keep this question in mind for Radjaï in Alert school. ;-)

Still, it is an "inertia" (if you don't like mass...) because it acts on
acceleration, while a stiffness would act on position. Stiffness and
mass are all "gains" if they appear in a servo-control loop.
In the sense you use it here, gain is a very vague concept. Strictly
speaking, the unit of gain is dB and it is a function of frequency...

Bruno

On 31 May 2011 08:41, Chareyre <email address hidden> wrote:

> Your question #159394 on Yade changed:
> https://answers.launchpad.net/yade/+question/159394
>
> Status: Open => Answered
>
> Chareyre proposed the following answer:
>
> > This equation is very similar to the servo-control we use with rigid
> > boundaries; with the difference that there we also have a gain parameter
> > (let me call it that way, for me it does not matter if it is split into
> two
> > parameters other than one). So, why here the formula is correct but with
> > rigid walls we need the gain parameter?
> For rigid walls (or periodic && !dynCell), it reads
> displacement=stress_error/stiffness in both cases (a bit more
> sophisticated in fact, but this is the basic idea).
> For dynCell, it reads acceleration=stress_error/mass (idem).
> I don't see what you mean by "we need a gain". Since a gain is any
> factor that enter a servo-control, there is basically a gain in all
> three cases.
>
Ok, you are right to say that there is a gain in all cases. I was making a
comparison between this control in the periodic

strain_rate=1/scene->dt*(goal[axis]-stress[axis])*cellArea[axis]/(stiff[axis]>0?stiff[axis]:1.);

and this control in triaxial test

previousTranslation[wall] = (1-wallDamping)*translation*normal[wall] +
0.8*previousTranslation[wall];

To me they are conceptually similar (leaving the 0.8 coeff aside), with the
difference that there is no wallDamping in periodic case. Maybe I am
misunderstanding...

>
> >
> > It is not good to have two methods if one of them does not work (at least
> to
> > my experience, but I would be happy to hear from others that I am wrong).
> Am
> > referring here to what you call "static" method (in PeriIso file). If you
> > ever try to use it, you will see that it does not properly work. This is
> why
> > my question about gain parameter.
>
> I think Vaclav and Jan used it successfully, with a different
> constitutive models than yours.
>
Vaclav did not use it and Jan wrote his own controller, but I let him add a
comment on all of this (could not you add your opinion, Jan?).

> I know Jan also tried to adjust the stiffness during iterations, instead
> of just predicting. So in principle it could work for many constitutive
> laws.
> I can't comment more personally, I never tried.
> > Mass is not mass, given its current unit of measure which can be derived
> > from the equation. As I see it, it is really a gain parameter (and I
> still
> > do not see the advantage to use 1/mass when this could simply be called
> > gain, unless you have a method which estimates the mass parameter from a
> > physical point of view, do you?).
> No, I don't have a clear derivation yet. I'm only guessing it needs to
> include lengths somewhere. It's an opened question for me (and if I read
> this book correctly, it's not explained either
> http://www.iste.co.uk/index.php?f=a&ACTION=View&id=384). I suggest you
> keep this question in mind for Radjaï in Alert school. ;-)
>
No worries, I have a list of questions ready to ask :-) Hopefully it will be
interesting.

>
> Still, it is an "inertia" (if you don't like mass...) because it acts on
> acceleration, while a stiffness would act on position. Stiffness and
> mass are all "gains" if they appear in a servo-control loop.
>
Definitely.

> In the sense you use it here, gain is a very vague concept. Strictly
> speaking, the unit of gain is dB and it is a function of frequency...
>
Maybe am not precise when I use this naming, but it is because I am used to
read about it in such a way.

Cheers.
Chiara

>
> Bruno
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/159394/+confirm?answer_id=2
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/159394
>
> You received this question notification because you asked the question.
>

> To me they are conceptually similar (leaving the 0.8 coeff aside), with the
> difference that there is no wallDamping in periodic case. Maybe I am
> misunderstanding...
That's right. In the periodic case, the equation imposes damping=0. If
it were me, I would have include (1-damping) factor in the equation.
If it is unstable, such factor could help relaxing the system (or in
other words: reduce gain).

> Vaclav did not use it and Jan wrote his own controller, but I let him add a
> comment on all of this (could not you add your opinion, Jan?).
Mmmh... Let's see.

Bruno

On 31 May 2011 14:41, Chareyre <email address hidden> wrote:

> Your question #159394 on Yade changed:
> https://answers.launchpad.net/yade/+question/159394
>
> Status: Open => Answered
>
> Chareyre proposed the following answer:
>
> > To me they are conceptually similar (leaving the 0.8 coeff aside), with
> the
> > difference that there is no wallDamping in periodic case. Maybe I am
> > misunderstanding...
> That's right. In the periodic case, the equation imposes damping=0. If
> it were me, I would have include (1-damping) factor in the equation.
> If it is unstable, such factor could help relaxing the system (or in
> other words: reduce gain).
>
Good, you come with me then... So we are missing a relaxation factor there,
I believe.
Bruno, I had a look at your equation, again

strain_rate+=dampFactor*scene->dt* ( goal[axis]-stress[axis] ) /mass;

and we are missing a length here to be consistent with the unit of measure.
Radjai, in a similar equation he proposes in the book you mention, is in
fact considering a length divided by the mass parameter. Why do not we use
it as well? We could take the size of the cell as a reference (Radjai's
suggestion).

Chiara

>
> > Vaclav did not use it and Jan wrote his own controller, but I let him add
> a
> > comment on all of this (could not you add your opinion, Jan?).
> Mmmh... Let's see.
>
> Bruno
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/159394/+confirm?answer_id=4
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/159394
>
> You received this question notification because you asked the question.
>

> strain_rate+=dampFactor*scene->dt* ( goal[axis]-stress[axis] ) /mass;
>
> and we are missing a length here to be consistent with the unit of measure.
> Radjai, in a similar equation he proposes in the book you mention, is in
> fact considering a length divided by the mass parameter. Why do not we use
> it as well? We could take the size of the cell as a reference (Radjai's
> suggestion).
>
So I missed it (which page?).

B.

On 31 May 2011 15:31, Chareyre <email address hidden> wrote:

> Your question #159394 on Yade changed:
> https://answers.launchpad.net/yade/+question/159394
>
> Status: Open => Answered
>
> Chareyre proposed the following answer:
>
> > strain_rate+=dampFactor*scene->dt* ( goal[axis]-stress[axis] ) /mass;
> >
> > and we are missing a length here to be consistent with the unit of
> measure.
> > Radjai, in a similar equation he proposes in the book you mention, is in
> > fact considering a length divided by the mass parameter. Why do not we
> use
> > it as well? We could take the size of the cell as a reference (Radjai's
> > suggestion).
> >
> So I missed it (which page?).
>
Page 168, hope I read well... Let me know what you think.

>
> B.
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/159394/+confirm?answer_id=6
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/159394
>
> You received this question notification because you asked the question.
>

>Vaclav did not use it and Jan wrote his own controller, but I let him add a
>comment on all of this (could not you add your opinion, Jan?).

Of course I could :-) but as I don't know PeriTriaxController and similar engines very well and having not so much time just now, could you please summarize your discussion here for me (if you still want my opinion)?

in the discussion there was mentioned using stiffness.. for stress/strain control of periodic cell, we tried to compute stiffness of whole cell according to published formulas, but the formulas itself were only approximation, so we found this method unsuitable. In current Peri3dController, the strain can be prescribed directly. For prescribed stress, the appropriate strain component is estimated so as the stress differs minimally from ideal value..

good luck in your work wishes
Jan

On 31 May 2011 17:35, honzik <email address hidden> wrote:

> Your question #159394 on Yade changed:
> https://answers.launchpad.net/yade/+question/159394
>
> Status: Open => Answered
>
> honzik proposed the following answer:
>
> >Vaclav did not use it and Jan wrote his own controller, but I let him add
> a
> >comment on all of this (could not you add your opinion, Jan?).
>
> Of course I could :-) but as I don't know PeriTriaxController and
> similar engines very well and having not so much time just now, could
> you please summarize your discussion here for me (if you still want my
> opinion)?
>
That is fine, Jan. It proves that I was right to say that the static method
was not used before (afaik). The line in question is

(in PeriTriaxController)
strain_rate=1/scene->dt*(goal[axis]-stress[axis])*cellArea[axis]/(stiff[axis]>0?stiff[axis]:1.);

To me, this needs to be adjusted. Am trying to do something with that now, I
hope that nobody will mind if I change something in the pertinent code (for
the static solution, at least).
Cheers.
Chiara

>
> in the discussion there was mentioned using stiffness.. for
> stress/strain control of periodic cell, we tried to compute stiffness of
> whole cell according to published formulas, but the formulas itself were
> only approximation, so we found this method unsuitable. In current
> Peri3dController, the strain can be prescribed directly. For prescribed
> stress, the appropriate strain component is estimated so as the stress
> differs minimally from ideal value..
>

>
> good luck in your work wishes
> Jan
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/159394/+confirm?answer_id=8
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/159394
>
> You received this question notification because you asked the question.
>

On 3 June 2011 09:10, Chareyre <email address hidden> wrote:

> Your question #159394 on Yade changed:
> https://answers.launchpad.net/yade/+question/159394
>
> Status: Open => Answered
>
> Chareyre proposed the following answer:
>
> >>> it as well? We could take the size of the cell as a reference (Radjai's
> >>> suggestion).
> >>>
> >> So I missed it (which page?).
> >>
> > Page 168, hope I read well... Let me know what you think.
> >
>
> Page 168 was written by my colleagues Gaël Combe and Jean-Noël Roux
> (while I was searching cell's mass in Radjai's PBC chapter) .
>
Oh sure, sorry if I confused the names.

> I got a chance to discuss the issue with them this week. Actually, when
> I asked which mass should be used, the first reply was "use a value that
> makes it work"... ;-)
> But still, this chapter suggest a definition that is dimensionally
> correct, you could try it.
> There is no need to change anything in the code, since it's up to you to
> give the "mass" value in your scripts. If you find that this definition
> is well suited, we can maybe hardcode it in the cpp later, so that
> "mass" will be totalMass/L by default (and include a mass factor to let
> users modify it if needed).
>
> Independently, I would keep the damping factor as it is now, since in
> perfectly elastic situations, eq. [6.7] would oscillate forever.
>
OK.

>
> In the end it would give:
>
> strain_rate[axis] += dampFactor*dt* ( goal[axis]-stress[axis] )
> *Length/(totalMass*massFactor); //massFactor=1 by default
>
I will try it. My conclusion is that probably it does not really matter how
you actually define your gain parameter, given that it works. However, it is
obviously good to have an idea (estimate) of what it could be and the above
line looks good to me. I will let you know. Thanks for your comments on
this.

Chiara

>
> Let us know if it works for you.
>
> Cheers.
>
> Bruno
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/159394/+confirm?answer_id=10
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/159394
>
> You received this question notification because you asked the question.
>

Thanks Chareyre, that solved my question.