Hi all,
I developed a "copy" of GlobalStiffnessTimeStepper which evaluate the time step
of a system of particle with Law2_ScGeom_ViscElPhys_Basic contact law (spring +
viscous damping).
The principle of the time step evaluation is :
We consider that the elastic and the viscous degrees of freedom are decoupled.
From that, we evaluate the critical time step of the system of particles
considering the elastic part only (dt_el) on one hand, and considering the
viscous part only (dt_visc) on the other hand. Then, we take the minimum of the
two timestep and multiply it by a safety factor to account for the possible
coupling between the two effect (dt = alpha*min(dt_el,dt_visc).
To do so, I took GlobalStiffnessTimeStepper (GSTS) which is already evaluating
the elastic critical time step, and implemented similarly a calculation for the
viscous critical time step of the system from the the viscosity matrix of each
particles. The changes are not so important between the version I implemented
and the last one (see the code below). I tested it in cases where the normal
viscous part is dominant and the new version gives a lower time step which
allows to keep the caculus stable, while in the old version everything explode,
so it seems to work.
I have some questions about that :
- Does the evaluation of the time step seems reasonable to you?
- What can I do to test it quantitavely ? (if I try a collision between two
beads, GSTS does not work..)
- To add it to the code, should I add it as a new "file"/"function" or should I
implement it as an option of GSTS ?
And about GSTS in general :
- why is GSTS not working when we consider a gravity deposition or a collision
?? (even when defaultDt is small enough)
Thanks,
Raphael
GlobalStiffnessTimeStepper.cpp
/*************************************************************************
* Copyright (C) 2006 by Bruno Chareyre *
* bruno.chare...@hmg.inpg.fr *
* *
* This program is free software; it is licensed under the terms of the *
* GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
#include"GlobalStiffnessTimeStepper.hpp"
#include<yade/pkg/dem/FrictPhys.hpp>
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/pkg/dem/DemXDofGeom.hpp>
#include<yade/core/Interaction.hpp>
#include<yade/core/Scene.hpp>
#include<yade/core/Clump.hpp>
#include<yade/pkg/dem/Shop.hpp>
#include<yade/pkg/dem/ViscoelasticPM.hpp>
CREATE_LOGGER(GlobalStiffnessTimeStepper);
YADE_PLUGIN((GlobalStiffnessTimeStepper));
GlobalStiffnessTimeStepper::~GlobalStiffnessTimeStepper() {}
void GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>&
body, Scene * ncb)
{
State* sdec=body->state.get();
Vector3r& stiffness= stiffnesses[body->getId()];
Vector3r& Rstiffness=Rstiffnesses[body->getId()];
Vector3r& viscosity= viscosities[body->getId()];
Vector3r& Rviscosity= Rviscosities[body->getId()];
if(body->isClump()) {// if clump, we sum stifnesses of all members
const shared_ptr<Clump>&
clump=YADE_PTR_CAST<Clump>(body->shape);
FOREACH(Clump::MemberMap::value_type& B, clump->members){
const shared_ptr<Body>& b = Body::byId(B.first,scene);
stiffness+=stiffnesses[b->getId()];
Rstiffness+=Rstiffnesses[b->getId()];
viscosity+=viscosities[b->getId()];
Rviscosity+=Rviscosities[b->getId()];
}
}
if(!sdec || stiffness==Vector3r::Zero()){
if (densityScaling) sdec->densityScaling =
min(1.0001*sdec->densityScaling,
timestepSafetyCoefficient*pow(defaultDt/targetDt,2.0));
return; // not possible to compute!
}
//Determine the elastic critical time step and the viscous one
separately (as if they were decoupled) and take the minimum time step with
a safety coefficient.
//Elastic
Real dtx, dty, dtz;
Real dt = max( max (stiffness.x(), stiffness.y()), stiffness.z() );
if (dt!=0) {
dt = sdec->mass/dt; computedSomething = true;}//dt = squared
eigenperiod of translational motion
else dt = Mathr::MAX_REAL;
if (Rstiffness.x()!=0) {
dtx = sdec->inertia.x()/Rstiffness.x(); computedSomething =
true;}//dtx = squared eigenperiod of rotational motion around x
else dtx = Mathr::MAX_REAL;
if (Rstiffness.y()!=0) {
dty = sdec->inertia.y()/Rstiffness.y(); computedSomething =
true;}
else dty = Mathr::MAX_REAL;
if (Rstiffness.z()!=0) {
dtz = sdec->inertia.z()/Rstiffness.z(); computedSomething =
true;}
else dtz = Mathr::MAX_REAL;
Real Rdt = std::min( std::min (dtx, dty), dtz );//Rdt = smallest
squared eigenperiod for elastic rotational motions
dt =
1.41044*timestepSafetyCoefficient*std::sqrt(std::min(dt,Rdt));//1.41044 =
sqrt(2)
//Viscous
Real dtx_visc, dty_visc, dtz_visc;
Real dt_visc = max(max(viscosity.x(), viscosity.y()), viscosity.z() );
if (dt_visc!=0) {
dt_visc = sdec->mass/dt_visc; computedSomething = true;}//dt =
eigenperiod of the viscous translational motion
else {dt_visc = Mathr::MAX_REAL;}
if (Rviscosity.x()!=0) {
dtx_visc = sdec->inertia.x()/Rviscosity.x(); computedSomething
= true;}//dtx = eigenperiod of viscous rotational motion around x
else dtx_visc = Mathr::MAX_REAL;
if (Rviscosity.y()!=0) {
dty_visc = sdec->inertia.y()/Rviscosity.y(); computedSomething
= true;}
else dty_visc = Mathr::MAX_REAL;
if (Rviscosity.z()!=0) {
dtz_visc = sdec->inertia.z()/Rviscosity.z(); computedSomething
= true;}
else dtz_visc = Mathr::MAX_REAL;
Real Rdt_visc = std::min( std::min (dtx_visc, dty_visc), dtz_visc
);//Rdt = smallest squared eigenperiod for viscous rotational motions
dt_visc = 2*timestepSafetyCoefficient*std::min(dt_visc,Rdt_visc);
//Take the minimum between the elastic and viscous critical time step.
dt = std::min(dt,dt_visc);
//if there is a target dt, then we apply density scaling on the body,
the inertia used in Newton will be mass*scaling, the weight is unmodified
if (densityScaling) {
sdec->densityScaling =
min(sdec->densityScaling,timestepSafetyCoefficient*pow(dt /targetDt,2.0));
newDt=targetDt;
}
//else we update dt normaly
else {newDt = std::min(dt,newDt);}
}
bool GlobalStiffnessTimeStepper::isActivated()
{
return (active && ((!computedOnce) || (scene->iter %
timeStepUpdateInterval == 0) || (scene->iter < (long int) 2) ));
}
void GlobalStiffnessTimeStepper::computeTimeStep(Scene* ncb)
{
// for some reason, this line is necessary to have correct functioning
(no idea _why_)
// see scripts/test/compare-identical.py, run with or without
active=active.
active=active;
if (defaultDt<0) defaultDt=
timestepSafetyCoefficient*Shop::PWaveTimeStep(Omega::instance().getScene());
computeStiffnesses(ncb);
shared_ptr<BodyContainer>& bodies = ncb->bodies;
newDt = Mathr::MAX_REAL;
computedSomething=false;
BodyContainer::iterator bi = bodies->begin();
BodyContainer::iterator biEnd = bodies->end();
for( ; bi!=biEnd ; ++bi ){
shared_ptr<Body> b = *bi;
if (b->isDynamic() && !b->isClumpMember())
findTimeStepFromBody(b, ncb);
}
if(densityScaling) (newDt=targetDt);
if(computedSomething || densityScaling){
previousDt = min ( min(newDt , maxDt), 1.05*previousDt );// at
maximum, dt will be multiplied by 1.05 in one iterration, this is to prevent
brutal switches from 0.000... to 1 in some computations
scene->dt=previousDt;
computedOnce = true;}
else if (!computedOnce) scene->dt=defaultDt;
LOG_INFO("computed timestep " << newDt <<
(scene->dt==newDt ? string(", applied") :
string(", BUT timestep is
")+lexical_cast<string>(scene->dt))<<".");
}
void GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
/* check size */
size_t size=stiffnesses.size();
if(size<rb->bodies->size()){
size=rb->bodies->size();
stiffnesses.resize(size); Rstiffnesses.resize(size);
viscosities.resize(size); Rviscosities.resize(size);
}
/* reset stored values */
memset(& stiffnesses[0],0,sizeof(Vector3r)*size);
memset(&Rstiffnesses[0],0,sizeof(Vector3r)*size);
memset(& viscosities[0],0,sizeof(Vector3r)*size);
memset(& Rviscosities[0],0,sizeof(Vector3r)*size);
FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
if(!contact->isReal()) continue;
GenericSpheresContact*
geom=YADE_CAST<GenericSpheresContact*>(contact->geom.get()); assert(geom);
ViscElPhys* phys = YADE_CAST<ViscElPhys*>(contact->phys.get());
assert(phys);
// all we need for getting stiffness
Vector3r& normal=geom->normal; Real& kn=phys->kn; Real&
ks=phys->ks; Real& radius1=geom->refR1; Real& radius2=geom->refR2;
Real& cn=phys->cn; Real& cs=phys->cs;
Real fn = (static_cast<NormShearPhys *>
(contact->phys.get()))->normalForce.squaredNorm();
if (fn==0) continue;//Is it a problem with some laws? I still
don't see why.
//Diagonal terms of the translational stiffness matrix
Vector3r diag_stiffness =
Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
diag_stiffness *= kn-ks;
diag_stiffness = diag_stiffness + Vector3r(1,1,1)*ks;
//diagonal terms of the rotational stiffness matrix
// Vector3r branch1 =
currentContactGeometry->normal*currentContactGeometry->radius1;
// Vector3r branch2 =
currentContactGeometry->normal*currentContactGeometry->radius2;
Vector3r diag_Rstiffness =
Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
std::pow(normal.x(),2)+std::pow(normal.z(),2),
std::pow(normal.x(),2)+std::pow(normal.y(),2));
diag_Rstiffness *= ks;
//Diagonal terms of the translational viscous matrix
Vector3r diag_viscosity =
Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
diag_viscosity *= cn-cs;
diag_viscosity = diag_viscosity + Vector3r(1,1,1)*cs;
//diagonal terms of the rotational viscous matrix
Vector3r diag_Rviscosity =
Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
std::pow(normal.x(),2)+std::pow(normal.z(),2),
std::pow(normal.x(),2)+std::pow(normal.y(),2));
diag_Rviscosity *= cs;
//NOTE : contact laws with moments would be handled correctly
by summing directly bending+twisting stiffness to diag_Rstiffness. The fact
that there is no problem currently with e.g. cohesiveFrict law is probably
because final computed dt is constrained by translational motion, not rotations.
stiffnesses [contact->getId1()]+=diag_stiffness;
Rstiffnesses[contact->getId1()]+=diag_Rstiffness*pow(radius1,2);
viscosities [contact->getId1()]+=diag_viscosity;
Rviscosities[contact->getId1()]+=diag_Rviscosity*pow(radius1,2);
stiffnesses [contact->getId2()]+=diag_stiffness;
Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
viscosities [contact->getId2()]+=diag_viscosity;
Rviscosities[contact->getId2()]+=diag_Rviscosity*pow(radius2,2);
}
}
GlobalStiffnessTimeStepper.hpp
/*************************************************************************
* Copyright (C) 2006 by Bruno Chareyre *
* bruno.chare...@hmg.inpg.fr *
* *
* This program is free software; it is licensed under the terms of the *
* GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
#pragma once
#include<yade/core/TimeStepper.hpp>
/*! \brief Compute the critical timestep of the leap-frog scheme based on
global stiffness of bodies
See usage details in TriaxialTest documentation (TriaxialTest is also a
good example of how to use this class)
*/
class Interaction;
class BodyContainer;
class Scene;
class GlobalStiffnessTimeStepper : public TimeStepper
{
private :
vector<Vector3r> stiffnesses;
vector<Vector3r> Rstiffnesses;
vector<Vector3r> viscosities;
vector<Vector3r> Rviscosities;
void computeStiffnesses(Scene*);
Real newDt;
bool computedSomething,
computedOnce;
void findTimeStepFromBody(const shared_ptr<Body>& body, Scene *
ncb);
public :
virtual ~GlobalStiffnessTimeStepper();
virtual void computeTimeStep(Scene*);
virtual bool isActivated();
YADE_CLASS_BASE_DOC_ATTRS_CTOR(
GlobalStiffnessTimeStepper,TimeStepper,"An engine
assigning the time-step as a fraction of the minimum eigen-period in the
problem. The derivation is detailed in the chapter on `DEM formulation
<https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_",
((Real,defaultDt,-1,,"used as the initial value of the
timestep (especially useful in the first steps when no contact exist). If
negative, it will be defined by :yref:`utils.PWaveTimeStep` *
:yref:`GlobalStiffnessTimeStepper::timestepSafetyCoefficient`"))
((Real,maxDt,Mathr::MAX_REAL,,"if positive, used as max
value of the timestep whatever the computed value"))
((Real,previousDt,1,,"last computed dt |yupdate|"))
((Real,timestepSafetyCoefficient,0.8,,"safety factor
between the minimum eigen-period and the final assigned dt (less than 1))"))
((bool,densityScaling,false,,"|yupdate| don't modify
this value if you don't plan to modify the scaling factor manually for some
bodies. In most cases, it is enough to set
:yref:`NewtonIntegrator::densityScaling` and let this one be adjusted
automatically."))
((Real,targetDt,1,,"if
:yref:`NewtonIntegrator::densityScaling` is active, this value will be used as
the simulation timestep and the scaling will use this value of dt as the
target value. The value of targetDt is arbitrary and should have no effect in
the result in general. However if some bodies have imposed velocities, for
instance, they will move more or less per each step depending on this value.")),
computedOnce=false;)
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(GlobalStiffnessTimeStepper);
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : yade-dev@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
