Skip to content

GitLab

Commit 708454b7 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

Cleaning scontact, removing relaxation and adding use_queue

parent 26646290
Pipeline #4941 passed with stage
in 26 seconds
Hide whitespace changes
Inline Side-by-side
python/scontact.py
View file @ 708454b7
......@@ -191,15 +191,12 @@ class ParticleProblem :
assert self._friction_enabled
self._lib.particleProblemSetBoundaryTangentVelocity(self._p,tag.encode(),c_double(v))
def set_friction_relaxation(self, relaxation):
""" Sets the relaxation parameter for the resolution of the frictional problem.
Only in 2D when friction is enabled.
Keyword arguments:
relaxation -- relaxation parameter, between 0 and 1
"""
assert self._friction_enabled
self._lib.particleProblemSetFrictionRelaxation(self._p,c_double(relaxation))
def set_use_queue(self, use_queue = True):
""" Sets the use of the queue in the resolution of the contacts."""
if use_queue:
self._lib.particleProblemSetUseQueue(self._p,c_int(1))
else:
self._lib.particleProblemSetUseQueue(self._p,c_int(0))
def set_friction_coefficient(self, mu, mat0="default",mat1="default") :
""" Sets the friction coefficient between two materials.
......@@ -222,7 +219,7 @@ class ParticleProblem :
tol -- optional argument defining the interpenetration tolerance to stop the NLGS iterations of the NSCD
"""
self._get_matrix("ExternalForces",self._dim)[:] = forces
return self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
def write_vtk(self, odir, i, t) :
"""Write output files for post-visualization."""
......
scontact/scontact.c
View file @ 708454b7
......@@ -117,8 +117,8 @@ static int particleInitContact(size_t id0, Particle *p0, double *x0, size_t id1,
c->a1 = p0->m / (p0->m + p1->m);
#ifdef FRICTION_ENABLED
c->ct = 0;
c->In0 = c->a0*5/(2*p0->r);
c->In1 = c->a1*5/(2*p1->r);
c->In0 = c->a0*2/(p0->r);
c->In1 = c->a1*2/(p1->r);
#endif
c->type = PARTICLE_PARTICLE;
return c->D < alert;
......@@ -181,7 +181,7 @@ static int diskInitContact(size_t id, const Disk *d, size_t particle, const Part
#ifdef FRICTION_ENABLED
c->ct = 0;
c->In0 = 0;
c->In1 = c->a1*5/(2*p->r);
c->In1 = c->a1*2/(p->r);
#endif
c->type = PARTICLE_DISK;
return c->D < alert;
......@@ -248,7 +248,7 @@ static int segmentInitContact(size_t id, const Segment *s, size_t particle, cons
#ifdef FRICTION_ENABLED
c->ct = 0;
c->In0 = 0;
c->In1 = c->a1*5/(2*p->r);
c->In1 = c->a1*2/(p->r);
#endif
c->type = PARTICLE_SEGMENT;
return c->D >= 0 && c->D < alert;
......@@ -336,7 +336,7 @@ static int triangleInitContact(size_t id, const Triangle *t, size_t particle, co
#ifdef FRICTION_ENABLED
c->ct = 0;
c->In0 = 0;
c->In1 = c->a1*5/(2*p->r);
c->In1 = c->a1*2/(p->r);
#endif
if (c->D < 0)
c->D = 0;
......@@ -362,7 +362,7 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
t[0] = c->n[1];
t[1] = -c->n[0];
double ct = 0 ;
double vt = //c->ct +
double vt =
c->ct*p->particles[c->o0].r*c->In0+
c->ct*p->particles[c->o1].r*c->In1+
(p->velocity[c->o0*2+0]-p->velocity[c->o1*2+0])*t[0]+
......@@ -370,18 +370,11 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
p->omega[c->o0]*p->particles[c->o0].r +
p->omega[c->o1]*p->particles[c->o1].r;
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o0],p->particleMaterial[c->o1]);
if(vt>(7./2.)*mu*(dp)){
ct = p->friction_relaxation*mu*(dp);
}
else if(vt<(-7./2.)*mu*(dp)){
ct = -p->friction_relaxation*mu*(dp);
}
else{
ct = p->friction_relaxation*(2./7.)*vt;
}
if (fabs(vt)> 3*mu*dp) vt *= 3*mu*dp/fabs(vt);
ct = vt/3;
if (dp == 0.)
ct = 0;
ct -= c->ct;
//coordAdd(&p->velocity[c->o0 * DIMENSION], -ct*c->a0, t);
//coordAdd(&p->velocity[c->o1 * DIMENSION], ct*c->a1, t);
p->omega[c->o0] -= c->In0*ct;
p->omega[c->o1] -= c->In1*ct;
c->ct += ct;
......@@ -392,7 +385,7 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
c->dv += dp;
if(fabs(dp) > tol/dt) return 0;
#ifdef FRICTION_ENABLED
if(fabs(ct/p->friction_relaxation) > tol/dt) return 0;
if(fabs(ct) > tol/dt) return 0;
#endif
return 1;
}
......@@ -415,7 +408,7 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d
double t[2];
t[0] = c->n[1];
t[1] = -c->n[0];
double vt = // c->ct +
double vt =
c->ct*p->particles[c->o1].r*c->In1+
p->velocity[c->o1 * DIMENSION]*t[0]+
p->velocity[c->o1 * DIMENSION+1]*t[1]+
......@@ -423,17 +416,11 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d
p->disks[c->o0].vt;
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->diskMaterial[c->o0]);
double ct = 0;
if(vt>(7./2.)*mu*dp){
ct = p->friction_relaxation*mu*dp;
}
else if(vt<(-7./2.)*mu*dp){
ct = -p->friction_relaxation*mu*dp;
}
else{
ct = p->friction_relaxation*(2./7.)*vt;
}
if (fabs(vt)> 3*mu*dp) vt *= 3*mu*dp/fabs(vt);
ct = vt/3;
if (dp == 0.)
ct = 0;
ct -= c->ct;
//coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*c->a1, t);
p->omega[c->o1] -= c->In1*ct;
c->ct += ct;
#endif
......@@ -442,7 +429,7 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d
c->dv += dp;
if(fabs(dp) > tol/dt) return 0;
#ifdef FRICTION_ENABLED
if(fabs(ct/p->friction_relaxation) > tol/dt) return 0;
if(fabs(ct) > tol/dt) return 0;
#endif
return 1;
}
......@@ -457,27 +444,23 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double dt
xn[i] = p->position[c->o1*DIMENSION + i] + r*c->n[i] + vlocfree[i]*c->D / vn;
if (!segmentProjectionIsInside(&p->segments[c->o0], xn))
dp = 0;
//if(c->o1 == 11)
//printf("contact %i %i , dp %g\n", c->o0, c->o1, dp);
#ifdef FRICTION_ENABLED
double t[2];
t[0] = c->n[1];
t[1] = -c->n[0];
double vt = // c->ct +
double vt =
c->ct*p->particles[c->o1].r*c->In1+
p->velocity[c->o1*DIMENSION]*t[0]+
p->velocity[c->o1*DIMENSION +1]*t[1]+
p->omega[c->o1]*p->particles[c->o1].r+
p->segments[c->o0].vt;
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->segmentMaterial[c->o0]);
printf("mu = %g\n", mu);
double ct = 0;
//if (fabs(vt)> 3.5*mu*dp) vt *= 3.5*mu*dp/fabs(vt);
ct = p->friction_relaxation*vt/3.5;
if (fabs(vt)> 3*mu*dp) vt *= 3*mu*dp/fabs(vt);
ct = vt/3;
if (dp == 0.)
ct = 0;
ct -= c->ct;
//coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*c->a1, t);
p->omega[c->o1] -= c->In1*ct;
c->ct += ct;
#endif
......@@ -485,12 +468,10 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double dt
coordAdd(&p->velocity[c->o1 * DIMENSION], -dp, c->n);
c->dv += dp;
#if FRICTION_ENABLED
printf("%g %g\n", t[0],t[1]);
printf("after solving %i %i : omega = %g, v = %g %g\n", c->o0, c->o1, p->omega[c->o1],p->velocity[c->o1 * DIMENSION+0],p->velocity[c->o1 * DIMENSION+1] );
#endif
if(fabs(dp) > tol/dt) return 0;
#if FRICTION_ENABLED
if(fabs(ct/p->friction_relaxation) > tol/dt) return 0;
if(fabs(ct) > tol/dt) return 0;
#endif
return 1;
}
......@@ -845,12 +826,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
coordAdd(&p->velocity[c->o1 * DIMENSION], cold->dv * c->a1, c->n);
c->dv = cold->dv;
#if FRICTION_ENABLED
double t[2];
t[0] = c->n[1];
t[1] = - c->n[0];
c->ct = cold->ct;
//coordAdd(&p->velocity[c->o0 * DIMENSION], -c->ct*c->a0, t);
//coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*c->a1, t);
p->omega[c->o0] -= c->In0*(c->ct);
p->omega[c->o1] -= c->In1*(c->ct);
#endif
......@@ -873,11 +849,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
coordAdd(&p->velocity[c->o1 * DIMENSION], -cold->dv, c->n);
c->dv = cold->dv;
#if FRICTION_ENABLED
double t[2];
t[0] = c->n[1];
t[1] = - c->n[0];
c->ct = cold->ct;
//coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t);
p->omega[c->o1] -= c->In1*(c->ct);
#endif
}
......@@ -897,11 +869,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
coordAdd(&p->velocity[c->o1 * DIMENSION], -cold->dv, c->n);
c->dv = cold->dv;
#if FRICTION_ENABLED
double t[2];
t[0] = c->n[1];
t[1] = -c->n[0];
c->ct = cold->ct;
//coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t);
p->omega[c->o1] -= c->In1*(c->ct);
#endif
}
......@@ -962,41 +930,19 @@ static void fifoFree(Fifo *f) {
free(f);
}
static int _particleProblemSolveContacts(ParticleProblem *p, double dt, double tol) {
static void _particleProblemSolveContacts(ParticleProblem *p, double dt, double tol) {
int converged = 0;
int iter = 0;
while (! converged) {
if (iter >= 0) {
#if FRICTION_ENABLED
printf("-----------------\n");
printf("----particles----\n");
for (int i = 0; i < vectorSize(p->particles); ++i) {
printf("%24.16e %24.16e -- %24.16e\n", p->position[i*2+0]+p->velocity[i*2+0]*dt,p->position[i*2+1]+p->velocity[i*2+1]*dt,p->omega[i]);
}
printf("----contacts----\n");
for (int i = 0; i < vectorSize(p->contacts); ++i) {
const Contact *c = &p->contacts[i];
printf("%3i %3i %3i : %24.16e %24.16e\n", c->type,c->o0,c->o1, c->dv, c->ct);
}
#endif
}
if(iter == 50) {
printf("blocked\n");
return 1;
}
iter++;
converged = 1;
for (int ic = vectorSize(p->contacts)-1; ic >= 0; --ic) {
Contact *c = &p->contacts[ic];
//if (iter%2 == 0)
//c = &p->contacts[vectorSize(p->contacts)-1-ic];
int conv;
switch (c->type) {
case PARTICLE_PARTICLE :
conv = contactParticleParticleSolve(c, p, dt,tol);
break;
case PARTICLE_DISK :
//conv = contactParticleDiskSolve(c, p, dt, tol);
conv = contactParticleDiskSolve(c, p, dt, tol);
break;
case PARTICLE_SEGMENT :
conv = contactParticleSegmentSolve(c, p, dt, tol);
......@@ -1012,7 +958,6 @@ static int _particleProblemSolveContacts(ParticleProblem *p, double dt, double t
}
}
}
return 0;
}
static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, double tol) {
......@@ -1104,13 +1049,13 @@ double particleProblemMaxDt(const ParticleProblem *p, double alert) {
return alert / q;
}
int particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit)
void particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit)
{
_particleProblemGenContacts(p, alert);
/*if (p->use_queue)
if (p->use_queue)
_particleProblemSolveContactsQueue(p,dt,tol);
else*/
return _particleProblemSolveContacts(p,dt,tol);
else
_particleProblemSolveContacts(p,dt,tol);
}
static void particleProblemApplyCtTranslation(ParticleProblem *p)
......@@ -1132,13 +1077,13 @@ static void particleProblemApplyCtTranslation(ParticleProblem *p)
}
int particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit)
void particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit)
{
int r = 0;
for (size_t j = 0; j < vectorSize(p->particles); ++j)
for (size_t i = 0; i < DIMENSION; ++i)
p->velocity[j * DIMENSION + i] += p->externalForces[j * DIMENSION + i] * dt / p->particles[j].m;
r = particleProblemSolve(p, alert, dt, tol, maxit);
particleProblemSolve(p, alert, dt, tol, maxit);
for (size_t i = 0; i < vectorSize(p->position); ++i){
p->position[i] += p->velocity[i] * dt;
}
......@@ -1160,8 +1105,6 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
}
#endif
particleProblemApplyCtTranslation(p);
printf("position : %g %g\n", p->position[0], p->position[1]);
return r;
}
size_t particleProblemAddBoundarySegmentTagId(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], int tag, int materialTag) {
......@@ -1402,7 +1345,6 @@ ParticleProblem *particleProblemNew() {
p->position = NULL;
p->velocity = NULL;
#ifdef FRICTION_ENABLED
p->friction_relaxation = 1.0;
p->omega = NULL;
p->mu = NULL;
#endif
......@@ -1449,9 +1391,6 @@ void particleProblemSetUseQueue(ParticleProblem *p, int use_queue) {
}
#ifdef FRICTION_ENABLED
void particleProblemSetFrictionRelaxation(ParticleProblem *p, double friction_relaxation) {
p->friction_relaxation = friction_relaxation;
}
double *particleProblemOmega(ParticleProblem *p) {return &p->omega[0];}
#endif
......
scontact/scontact.h
View file @ 708454b7
......@@ -31,8 +31,8 @@ ParticleProblem *particleProblemNew();
void particleProblemDelete(ParticleProblem *p);
void particleProblemLoad(ParticleProblem *p, const char *filename);
void particleProblemWrite(const ParticleProblem *p, const char *filename);
int particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit);
int particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit);
void particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit);
void particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit);
double particleProblemMaxDt(const ParticleProblem *p, double alert);
void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r, double m, const char *material);
size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x0[DIMENSION], double r, const char *tag, const char *material);
......@@ -68,7 +68,6 @@ void ParticleToMesh(size_t nParticles, double *particles, int nElements, double
#ifdef FRICTION_ENABLED
void particleProblemSetFrictionCoefficient(ParticleProblem *p, double mu, const char *mat0, const char *mat1);
double *particleProblemOmega(ParticleProblem *p);
void particleProblemSetFrictionRelaxation(ParticleProblem *p, double friction_relaxation);
#endif
int particleProblemNTag(const ParticleProblem *p);
......
testcases/mill-2d/mill.py
View file @ 708454b7
......@@ -64,15 +64,15 @@ ii = 0
#physical parameters
g = -9.81 # gravity
rhop = 2800 # grains density
tEnd = 1000 # final time
tEnd = 5 # final time
#numerical parameters
dt = 1e-3 # time step
#geometry parameters
rout = 0.11 # outer radius
r = 1.7e-2 # grains radius
Fr = 6 #Froude number
r = 1.8e-3 # grains radius
Fr = 0.5 #Froude number
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
v = rout*np.sqrt(Fr*(-g)/rout)
print('v = %s\n'%(v))
......@@ -92,9 +92,7 @@ outf = 1 # number of iterations between output fi
p.set_friction_coefficient(1.5,"sand","rock")
p.set_friction_coefficient(1.5,"sand","sand")
p.set_friction_coefficient(1.5,"rock","rock")
p.set_friction_relaxation(1.)
p.set_tangent_boundary_velocity("Outer",v)
print(p.segments()[:,:4])
#Computation loop
forces = np.zeros_like(p.velocity())
k=0
......@@ -112,7 +110,7 @@ while t < tEnd :
k = k+1
for i in range(nsub) :
tol = 1e-6
r = p.iterate(dt/nsub, forces, tol)
p.iterate(dt/nsub, forces, tol)
if r==1 :
ioutput = int(ii/outf) + 2
p.write_vtk(outputdir, ioutput, t)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment