Velocity gradient

I will have to think twice...
Is it related to the problem that Vaclav mentioned for homo=3?

Sorry, you lost me.

line 120 in Newton.cpp, I see:
scene->forces.addTorque(id,m);

line 178:
} YADE_PARALLEL_FOREACH_BODY_END();

I just synced with bzr.

Sorry! I got the wrong numbers. So here they are.

Lines 135:
Vector3r fluctVel=isPeriodic?scene->cell->bodyFluctuationVel(b->state->pos,b
->state->vel):state->vel;
and line 193:
if (homoDeform==Cell::HOMO_VEL_2ND) state->vel+=prevVelGrad*state->vel*dt;

It is not related to the problem of homoDeform, that was a mistake in the
code not in the formulation. So what do you think?
Chiara

On 31 August 2011 17:55, Chareyre <email address hidden>wrote:

> Your question #169689 on Yade changed:
> https://answers.launchpad.net/yade/+question/169689
>
> Chareyre proposed the following answer:
> Sorry, you lost me.
>
> line 120 in Newton.cpp, I see:
> scene->forces.addTorque(id,m);
>
> line 178:
> } YADE_PARALLEL_FOREACH_BODY_END();
>
> I just synced with bzr.
>
> --
> 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/169689/+confirm?answer_id=1
>
> 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/169689
>
> You received this question notification because you asked the question.
>

1) I agree to use prevVelGrad.

2) Indeed it's using v(t+dt/2) (more or less since it's not fully
updated). I think it should be using the value v(t-dt/2). How would it
read at time "t"?
I imagine this cycle where the body moves only due to cell deformation
(no contact forces), with velGrad switching between -0.5 and +0.5 at
each cycle, and assuming dt=1:

* position=1, velgrad=-0.5 : position goes to 0.5 after the first
integration, as velocity is -0.5
* next step, velgrad=+0.5 : velocity is adjusted at line 193:
v=v+0.25=-0.25 (correct, it corresponds to the velocity when the
particle arrives at position 0.5 in the old velocity field); then it's
accelerated the other way by dVelGrad=+2, v=0.25, pos=0.75.

The result after 2 steps looks correct (the position is not back to 1,
this is expected as the final size of the period itself is really 0.75,
not 1 (1*0.5*1.5)).
I commited a fix for 2), tell me if you agree (just moving a few lines
and updating acceleration instead of velocity).

On 1 September 2011 11:30, Chareyre <email address hidden>wrote:

> Your question #169689 on Yade changed:
> https://answers.launchpad.net/yade/+question/169689
>
> Status: Open => Answered
>
> Chareyre proposed the following answer:
> 1) I agree to use prevVelGrad.
>
Ok, let's do it then. We should simple add the velGrad argument to that
function and then input the prevVelGrad value when it is called in the
integrator.

>
> 2) Indeed it's using v(t+dt/2) (more or less since it's not fully
> updated). I think it should be using the value v(t-dt/2). How would it
> read at time "t"?
>
Well the correction to the velocity due to the homogeneous deformation is
(was, before last changes) made up of two terms.
First term (this is at time t, so it is fine): Vector3r dVel=dVelGrad*state
->pos;
Second term (it should have been at time t but it is/was not): state->vel+=
scene->cell->prevVelGrad*state->vel*dt;

The second term could be expressed like
(prevVelGrad+velGrad)/2 * (v(t+dt/2)-a(t)*dt/2) * dt

I shall look at what you propose and I let you know. You also let me know
what you think about it.

I imagine this cycle where the body moves only due to cell deformation
> (no contact forces), with velGrad switching between -0.5 and +0.5 at
> each cycle, and assuming dt=1:
>
> * position=1, velgrad=-0.5 : position goes to 0.5 after the first
> integration, as velocity is -0.5
> * next step, velgrad=+0.5 : velocity is adjusted at line 193:
> v=v+0.25=-0.25 (correct, it corresponds to the velocity when the
> particle arrives at position 0.5 in the old velocity field); then it's
> accelerated the other way by dVelGrad=+2, v=0.25, pos=0.75.
>
> The result after 2 steps looks correct (the position is not back to 1,
> this is expected as the final size of the period itself is really 0.75,
> not 1 (1*0.5*1.5)).
> I commited a fix for 2), tell me if you agree (just moving a few lines
> and updating acceleration instead of velocity).
>
> --
> 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/169689/+confirm?answer_id=3
>
> 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/169689
>
> You received this question notification because you asked the question.
>

For 1), you can also use cell->prevVelGrad in the function and
everywhere in Newton (and update it). This attribute is currently useless.

For 2), I think the two corrections should really be at different times.
The convective term is a correction of velocity accounting for the
displacement that occured between t-dt and t. This displacement is
exactly v(t-dt/2)*dt.
The term with dVelGrad reflects an "impulse" (or acceleration) occuring
at time t, yes, but it actually defines the velocity between t and t+dt
(what we call v(t+dt/2).
In other words, we correct the old, then we define the new.

Computing an average velGrad for time t would not make much sense since
velocities are never defined for integer steps.
If you apply this idea in the cycle I described, if I'm not wrong, it
would result in:

* first step: v=v-0.25*0=0 (convective); v=dVelGrad*pos=-0.5 (like before)
* second step: v=v+0*meanVel=-0.5 (no correction); v=v+dVelGrad*pos=0

As you see the particle is stuck at second step. The cell is deforming
but the particle doesn't move.
I'm not really sure I'm applying your suggestion correctly because in
fact this equation looks recursive to me:

v+=(prevVelGrad+velGrad)/2 * (v(t+dt/2)-a(t)*dt/2) * dt

It means:

v(t+dt/2)=v(t)+(prevVelGrad+velGrad)/2 * (v(t+dt/2)-a(t)*dt/2) * dt

How can we use v(t+dt/2) to define v(t+dt/2)?
Same remark for a(t), which is unknown before we apply this correction.

On 1 September 2011 17:00, Chareyre <email address hidden>wrote:

> Your question #169689 on Yade changed:
> https://answers.launchpad.net/yade/+question/169689
>
> Status: Open => Answered
>
> Chareyre proposed the following answer:
> For 1), you can also use cell->prevVelGrad in the function and
> everywhere in Newton (and update it). This attribute is currently useless.
>
> For 2), I think the two corrections should really be at different times.
> The convective term is a correction of velocity accounting for the
> displacement that occured between t-dt and t. This displacement is
> exactly v(t-dt/2)*dt.
> The term with dVelGrad reflects an "impulse" (or acceleration) occuring
> at time t, yes, but it actually defines the velocity between t and t+dt
> (what we call v(t+dt/2).
> In other words, we correct the old, then we define the new.
>
> Computing an average velGrad for time t would not make much sense since
> velocities are never defined for integer steps.
>
In fact I was considering to compute the two corrections to the velocity
both at time t. As you say, we use dVelGrad to update velocity at mid-step.
For me it was correct to consider _both_ corrections at time t, being the
value we update at time (t+dt/2). Anyway there is the problem below where I
would use v(t+dt/2) to update v(t+dt/2). In fact the latter is not really
v(t+dt/2) since this is not yet updated.

I can only partially see what you are doing by splitting these two
contributions, one acting on the acceleration and the other one on the
velocity. Reading the chapter by Radjai in the Discrete Element book about
PBCs, I can see that they apply ONE single corrective term to the
acceleration (please see eq. 7.26 at pag. 190). However there is a
difference between their formulation and ours since they distinguish the
fluctuating velocity of the fluid particle to the velocity field applied
whereas we do not do that if not when we compute kinetic energy.

Can we somehow be consistent with their formulation which seems good to me?

Chiara

> If you apply this idea in the cycle I described, if I'm not wrong, it
> would result in:
>
> * first step: v=v-0.25*0=0 (convective); v=dVelGrad*pos=-0.5 (like before)
> * second step: v=v+0*meanVel=-0.5 (no correction); v=v+dVelGrad*pos=0
>
> As you see the particle is stuck at second step. The cell is deforming
> but the particle doesn't move.
> I'm not really sure I'm applying your suggestion correctly because in
> fact this equation looks recursive to me:
>
> v+=(prevVelGrad+velGrad)/2 * (v(t+dt/2)-a(t)*dt/2) * dt
>
>
> It means:
>
> v(t+dt/2)=v(t)+(prevVelGrad+velGrad)/2 * (v(t+dt/2)-a(t)*dt/2) * dt
>
>
> How can we use v(t+dt/2) to define v(t+dt/2)?
> Same remark for a(t), which is unknown before we apply this correction.
>
> --
> 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/169689/+confirm?answer_id=5
>
> 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/169689
>
> You received this question notification because you asked the question.
>

For me eq. 7.26 is very close to:

Vector3r dVel=dVelGrad*state->pos;
state->vel+=dVel;

dVel is an impulse at time t right?

Now the other correction we have is either hidden in their equations,
due to the difference between fluctuational and total velocity, or due a
mistake in a derivation on one side.
I wouldn't bet my life on that but when I checked I found a mistake in
their derivation. I'll try and find it again and will let you know.

> Now the other correction we have is either hidden in their equations,
> due to the difference between fluctuational and total velocity, or due a
> mistake in a derivation on one side.
The first is true. In Radjai, the contact laws must compute the absolute
velocity = fluctuation+velgrad*pos (1).
In Yade we assign the absolute velocity in Newton, so that in contact
laws we don't care summing two terms.
In turn we have to add velG*pos and remove it at the next time-step.

You can see that by rewriting the operations differently. We have now:

v+=VG(t-1/2)*v(t-1/2)*dt
v+=(VG(t+1/2)-VG(t-1/2))*pos(t)

Due the fact that v(t-1/2)*dt=pos(t)-pos(t-1), it is equivalent to:

v+= -VG(t-1/2)*pos(t-1) + VG(t+1/2)*pos(t)

i.e. we remove the affine velocity of the previous step from the
absolute velocity (it leaves the fluctuation only) and we add the affine
velocity of the next step.
We are exactly equivalent I think.

When you say Radjai applies only one correction, it is not true. He has
only one correction in Newton, but he has another one in contact laws.
We have two in Newton but none in contact laws.
Since one body's velocity will be used for many contacts, it seems more
optimal to compute it once in Newton. I don't think we have to change
anything.