Commit 4fdfb5b1 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

possible solution for convergence

parent e0001b5b
Pipeline #6496 passed with stage
in 2 minutes and 9 seconds
......@@ -233,7 +233,7 @@ class ParticleProblem :
assert self._friction_enabled
self._lib.particleProblemSetFrictionCoefficient(self._p, c_double(mu), mat0.encode(),mat1.encode())
def iterate(self, dt, forces,tol=1e-8) :
def iterate(self, dt, forces,tol=1e-8,force_motion=0) :
"""Compute iteratively the set of velocities that respect the non-interpenetration constrain.
Keyword arguments:
......@@ -243,8 +243,10 @@ class ParticleProblem :
"""
vn = self.velocity().copy()
self._get_matrix("ExternalForces",self._dim)[:] = forces
self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
self._lib.particleProblemIterate.restype = c_int
retour = self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1),c_int(force_motion))
self._fc = (self.velocity()-vn)*self.mass()/dt - forces
return retour
def write_vtk(self, odir, i, t) :
"""Write output files for post-visualization."""
......
import numpy as np
def _advance_particles(particles, f, dt, min_nsub,contact_tol) :
def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0) :
# Compute free velocities of the grains
if f is None:
f = np.zeros_like(particles.velocity())
# Compute free velocities of the grains
vn = particles.velocity() + f*dt / particles.mass()
if particles.velocity().shape[1] == 2:
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
else:
if particles.velocity().shape[1] - 2:
vmax = np.max(np.sqrt(vn[:, 0]**2+ vn[:, 1]**2 + vn[:,2]**2))
else:
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
# Estimation of the solid sub-time steps
nsub = max(min_nsub, int(np.ceil((vmax * dt * 8)/min(particles.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(particles.r()))
# For each sub-time step solve contacts (do NLGS iterations) and move the grains
#Ignore the convergence and move the grains
if iteration == 10:
for i in range(nsub) :
particles.iterate(dt/nsub, f, tol=contact_tol,force_motion=1)
return
# For each sub-time step solve contacts (do NLGS iterations) and move the grains
for i in range(nsub) :
particles.iterate(dt/nsub, f, tol=contact_tol)
if not particles.iterate(dt/nsub, f, tol=contact_tol):
print("Splitting time-step to level %d..."%(iteration+1))
_advance_particles(particles,f,dt/nsub,2,contact_tol,iteration+1)
print("Convergence reached for subinvervals of level %d"%(iteration+1))
def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5) :
"""Predictor-corrector scheme to solve fluid and grains.
......
......@@ -332,22 +332,9 @@ static void _cross (const double *a, const double *b, double *c) {
}
static void build_basis(const double *a, double *b, double *c){
double max = fmax(fmax(fabs(a[0]), fabs(a[1])), fabs(a[2]));
if(1){
b[0] = -a[2];
b[1] = a[0];
b[2] = -a[1];
}
else if (fabs(a[1])==max){
b[0] = 0;
b[1] = a[2];
b[2] = -a[1];
}
else if(fabs(a[2])==max){
b[0] = -a[2];
b[1] = 0;
b[2] = a[1];
}
b[0] = -a[2];
b[1] = a[0];
b[2] = -a[1];
_cross(a, b, c);
}
......@@ -1371,14 +1358,14 @@ static void fifoFree(Fifo *f) {
free(f);
}
static void _particleProblemSolveContacts(ParticleProblem *p, double dt, double tol) {
static int _particleProblemSolveContacts(ParticleProblem *p, double dt, double tol) {
int converged = 0;
#if FRICTION_ENABLED
int iter = 0;
#endif
while (! converged) {
#if FRICTION_ENABLED
if(iter>100*vectorSize(p->particles)) break;
if(iter>100*vectorSize(p->particles)) return 0;
#endif
converged = 1;
for (int ic = vectorSize(p->contacts)-1; ic >= 0; --ic) {
......@@ -1408,9 +1395,10 @@ static void _particleProblemSolveContacts(ParticleProblem *p, double dt, double
iter++;
#endif
}
return 1;
}
static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, double tol) {
static int _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, double tol) {
size_t *nContactByParticle = malloc(sizeof(size_t)*(vectorSize(p->particles)));
for (size_t i = 0; i < vectorSize(p->particles); ++i) {
nContactByParticle[i] = 0;
......@@ -1451,28 +1439,7 @@ static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, do
#if FRICTION_ENABLED
if(iter>2000*vectorSize(p->contacts)*(DIMENSION)){
printf("Warning : Convergence not reached\n");
for (int ik = vectorSize(p->contacts)-1; ik >= 0; --ik) {
Contact *d = &p->contacts[ik];
int conv;
switch (d->type) {
case PARTICLE_PARTICLE :
conv = contactParticleParticleSolve(d, p, dt,tol);
break;
case PARTICLE_DISK :
conv = contactParticleDiskSolve(d, p, dt, tol);
break;
case PARTICLE_SEGMENT :
conv = contactParticleSegmentSolve(d, p, dt, tol);
break;
#if DIMENSION == 3
case PARTICLE_TRIANGLE :
conv = contactParticleTriangleSolve(d, p, dt, tol);
break;
#endif
}
}
break;
return 0;
}
#endif
......@@ -1514,7 +1481,7 @@ static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, do
#endif
activeContact[ic] = 0;
}
return 1;
fifoFree(queue);
free(activeContact);
free(contactByParticleP);
......@@ -1534,21 +1501,24 @@ double particleProblemMaxDt(const ParticleProblem *p, double alert) {
return alert / q;
}
void particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit)
int particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit)
{
_particleProblemGenContacts(p, alert);
if (p->use_queue)
_particleProblemSolveContactsQueue(p,dt,tol);
return _particleProblemSolveContactsQueue(p,dt,tol);
else
_particleProblemSolveContacts(p,dt,tol);
return _particleProblemSolveContacts(p,dt,tol);
}
void particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit)
int particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit,int force_motion)
{
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;
particleProblemSolve(p, alert, dt, tol, maxit);
if(!particleProblemSolve(p, alert, dt, tol, maxit) && !force_motion){
return 0;
}
for (size_t i = 0; i < vectorSize(p->position); ++i){
p->position[i] += p->velocity[i] * dt;
}
......@@ -1569,79 +1539,7 @@ void particleProblemIterate(ParticleProblem *p, double alert, double dt, double
}
}
#endif
if(p->get_momentum ==1){
for(size_t k = 0; k<DIMENSION*vectorSize(p->tagname);k++)
p->momentum[k] = 0.;
//printf("coucou\n");
for(size_t i = 0; i<vectorSize(p->contacts); ++i){
if(p->contacts[i].type == PARTICLE_DISK ){
double n[2];
n[0] = p->contacts[i].n[0];
n[1] = p->contacts[i].n[1];
#if FRICTION_ENABLED && DIMENSION==2
double t[2];
t[0] = -n[1];
t[1] = n[0];
#elif FRICTION_ENABLED && DIMENSION==3
double s[3], t[3];
t[0] = -n[1];
t[1] = n[0];
t[2] = n[2];
_cross(n, t, s);
#endif
for(int j = 0; j<DIMENSION;j++){
p->momentum[DIMENSION*p->diskTag[p->contacts[i].o0] +j] += n[j]*p->particles[p->contacts[i].o1].m*p->contacts[i].dv;
#if FRICTION_ENABLED && ROTATION_ENABLED
p->momentum[DIMENSION*p->diskTag[p->contacts[i].o0] +j] += t[j]*3*p->particles[p->contacts[i].o1].m*p->contacts[i].ct;
#if DIMENSION==3
p->momentum[DIMENSION*p->diskTag[p->contacts[i].o0] +j] += s[j]*3*p->particles[p->contacts[i].o1].m*p->contacts[i].cs;
#endif
#elif FRICTION_ENABLED
p->momentum[DIMENSION*p->diskTag[p->contacts[i].o0] +j] += t[j]*p->particles[p->contacts[i].o1].m*p->contacts[i].ct;
#if DIMENSION==3
p->momentum[DIMENSION*p->diskTag[p->contacts[i].o0] +j] += s[j]*p->particles[p->contacts[i].o1].m*p->contacts[i].cs;
#endif
#endif
}
}
else if (p->contacts[i].type == PARTICLE_SEGMENT){
double n[2];
n[0] = p->contacts[i].n[0];
n[1] = p->contacts[i].n[1];
#if FRICTION_ENABLED && DIMENSION==2
double t[2];
t[0] = -n[1];
t[1] = n[0];
#elif FRICTION_ENABLED && DIMENSION==3
double s[3], t[3];
t[0] = -n[1];
t[1] = n[0];
t[2] = n[2];
_cross(n, t, s);
#endif
for(int j = 0; j<DIMENSION;j++){
p->momentum[DIMENSION*p->segmentTag[p->contacts[i].o0] +j] += n[j]*p->particles[p->contacts[i].o1].m*p->contacts[i].dv;
#if FRICTION_ENABLED && ROTATION_ENABLED
p->momentum[DIMENSION*p->segmentTag[p->contacts[i].o0] +j] += t[j]*3*p->particles[p->contacts[i].o1].m*p->contacts[i].ct;
#if DIMENSION==3
p->momentum[DIMENSION*p->segmentTag[p->contacts[i].o0] +j] += s[j]*3*p->particles[p->contacts[i].o1].m*p->contacts[i].cs;
#endif
#elif FRICTION_ENABLED
p->momentum[DIMENSION*p->segmentTag[p->contacts[i].o0] +j] += t[j]*p->particles[p->contacts[i].o1].m*p->contacts[i].ct;
#if DIMENSION==3
p->momentum[DIMENSION*p->segmentTag[p->contacts[i].o0] +j] += s[j]*p->particles[p->contacts[i].o1].m*p->contacts[i].cs;
#endif
#endif
}
}
}
}
return 1;
}
size_t particleProblemAddBoundarySegmentTagId(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], int tag, int materialTag) {
......
......@@ -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);
void particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit);
void particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit);
int particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit,int force_motion);
int 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);
......
Supports Markdown
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