yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #09912
Modified GlobalStiffnessTimeStepper for visco-elastic contact law
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.chareyre@xxxxxxxxxxx *
* *
* 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.chareyre@xxxxxxxxxxx *
* *
* 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);
Follow ups
