Commit 26646290 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

wip decoupled

parent 70923058
Pipeline #4937 passed with stage
in 26 seconds
...@@ -222,7 +222,7 @@ class ParticleProblem : ...@@ -222,7 +222,7 @@ class ParticleProblem :
tol -- optional argument defining the interpenetration tolerance to stop the NLGS iterations of the NSCD tol -- optional argument defining the interpenetration tolerance to stop the NLGS iterations of the NSCD
""" """
self._get_matrix("ExternalForces",self._dim)[:] = forces 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)) return 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) : def write_vtk(self, odir, i, t) :
"""Write output files for post-visualization.""" """Write output files for post-visualization."""
......
...@@ -362,7 +362,7 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d ...@@ -362,7 +362,7 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
t[0] = c->n[1]; t[0] = c->n[1];
t[1] = -c->n[0]; t[1] = -c->n[0];
double ct = 0 ; double ct = 0 ;
double vt = c->ct + double vt = //c->ct +
c->ct*p->particles[c->o0].r*c->In0+ c->ct*p->particles[c->o0].r*c->In0+
c->ct*p->particles[c->o1].r*c->In1+ c->ct*p->particles[c->o1].r*c->In1+
(p->velocity[c->o0*2+0]-p->velocity[c->o1*2+0])*t[0]+ (p->velocity[c->o0*2+0]-p->velocity[c->o1*2+0])*t[0]+
...@@ -380,8 +380,8 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d ...@@ -380,8 +380,8 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
ct = p->friction_relaxation*(2./7.)*vt; ct = p->friction_relaxation*(2./7.)*vt;
} }
ct -= c->ct; ct -= c->ct;
coordAdd(&p->velocity[c->o0 * DIMENSION], -ct*c->a0, t); //coordAdd(&p->velocity[c->o0 * DIMENSION], -ct*c->a0, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], ct*c->a1, t); //coordAdd(&p->velocity[c->o1 * DIMENSION], ct*c->a1, t);
p->omega[c->o0] -= c->In0*ct; p->omega[c->o0] -= c->In0*ct;
p->omega[c->o1] -= c->In1*ct; p->omega[c->o1] -= c->In1*ct;
c->ct += ct; c->ct += ct;
...@@ -415,7 +415,7 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d ...@@ -415,7 +415,7 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d
double t[2]; double t[2];
t[0] = c->n[1]; t[0] = c->n[1];
t[1] = -c->n[0]; t[1] = -c->n[0];
double vt = c->ct + double vt = // c->ct +
c->ct*p->particles[c->o1].r*c->In1+ c->ct*p->particles[c->o1].r*c->In1+
p->velocity[c->o1 * DIMENSION]*t[0]+ p->velocity[c->o1 * DIMENSION]*t[0]+
p->velocity[c->o1 * DIMENSION+1]*t[1]+ p->velocity[c->o1 * DIMENSION+1]*t[1]+
...@@ -433,7 +433,7 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d ...@@ -433,7 +433,7 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d
ct = p->friction_relaxation*(2./7.)*vt; ct = p->friction_relaxation*(2./7.)*vt;
} }
ct -= c->ct; ct -= c->ct;
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*c->a1, t); //coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*c->a1, t);
p->omega[c->o1] -= c->In1*ct; p->omega[c->o1] -= c->In1*ct;
c->ct += ct; c->ct += ct;
#endif #endif
...@@ -457,37 +457,39 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double dt ...@@ -457,37 +457,39 @@ 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; xn[i] = p->position[c->o1*DIMENSION + i] + r*c->n[i] + vlocfree[i]*c->D / vn;
if (!segmentProjectionIsInside(&p->segments[c->o0], xn)) if (!segmentProjectionIsInside(&p->segments[c->o0], xn))
dp = 0; dp = 0;
//if(c->o1 == 11)
//printf("contact %i %i , dp %g\n", c->o0, c->o1, dp);
#ifdef FRICTION_ENABLED #ifdef FRICTION_ENABLED
double t[2]; double t[2];
t[0] = c->n[1]; t[0] = c->n[1];
t[1] = -c->n[0]; t[1] = -c->n[0];
double vt = c->ct + double vt = // c->ct +
c->ct*p->particles[c->o1].r*c->In1+ c->ct*p->particles[c->o1].r*c->In1+
p->velocity[c->o1*DIMENSION]*t[0]+ p->velocity[c->o1*DIMENSION]*t[0]+
p->velocity[c->o1*DIMENSION +1]*t[1]+ p->velocity[c->o1*DIMENSION +1]*t[1]+
p->omega[c->o1]*p->particles[c->o1].r+ p->omega[c->o1]*p->particles[c->o1].r+
p->segments[c->o0].vt; p->segments[c->o0].vt;
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->segmentMaterial[c->o0]); const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->segmentMaterial[c->o0]);
printf("mu = %g\n", mu);
double ct = 0; double ct = 0;
if(vt>(7./2.)*mu*dp){ //if (fabs(vt)> 3.5*mu*dp) vt *= 3.5*mu*dp/fabs(vt);
ct = p->friction_relaxation*mu*dp; ct = p->friction_relaxation*vt/3.5;
} if (dp == 0.)
else if(vt<(-7./2.)*mu*dp){ ct = 0;
ct = -p->friction_relaxation*mu*dp;
}
else{
ct = p->friction_relaxation*(2./7.)*vt;;
}
ct -= c->ct; ct -= c->ct;
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*c->a1, t); //coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*c->a1, t);
p->omega[c->o1] -= c->In1*ct; p->omega[c->o1] -= c->In1*ct;
c->ct += ct; c->ct += ct;
#endif #endif
dp -= c->dv; dp -= c->dv;
coordAdd(&p->velocity[c->o1 * DIMENSION], -dp, c->n); coordAdd(&p->velocity[c->o1 * DIMENSION], -dp, c->n);
c->dv += dp; 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(fabs(dp) > tol/dt) return 0;
#ifdef FRICTION_ENABLED #if FRICTION_ENABLED
if(fabs(ct/p->friction_relaxation) > tol/dt) return 0; if(fabs(ct/p->friction_relaxation) > tol/dt) return 0;
#endif #endif
return 1; return 1;
...@@ -847,8 +849,8 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert) ...@@ -847,8 +849,8 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
t[0] = c->n[1]; t[0] = c->n[1];
t[1] = - c->n[0]; t[1] = - c->n[0];
c->ct = cold->ct; c->ct = cold->ct;
coordAdd(&p->velocity[c->o0 * DIMENSION], -c->ct*c->a0, t); //coordAdd(&p->velocity[c->o0 * DIMENSION], -c->ct*c->a0, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*c->a1, t); //coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*c->a1, t);
p->omega[c->o0] -= c->In0*(c->ct); p->omega[c->o0] -= c->In0*(c->ct);
p->omega[c->o1] -= c->In1*(c->ct); p->omega[c->o1] -= c->In1*(c->ct);
#endif #endif
...@@ -875,7 +877,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert) ...@@ -875,7 +877,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
t[0] = c->n[1]; t[0] = c->n[1];
t[1] = - c->n[0]; t[1] = - c->n[0];
c->ct = cold->ct; c->ct = cold->ct;
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t); //coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t);
p->omega[c->o1] -= c->In1*(c->ct); p->omega[c->o1] -= c->In1*(c->ct);
#endif #endif
} }
...@@ -899,7 +901,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert) ...@@ -899,7 +901,7 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
t[0] = c->n[1]; t[0] = c->n[1];
t[1] = -c->n[0]; t[1] = -c->n[0];
c->ct = cold->ct; c->ct = cold->ct;
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t); //coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t);
p->omega[c->o1] -= c->In1*(c->ct); p->omega[c->o1] -= c->In1*(c->ct);
#endif #endif
} }
...@@ -960,19 +962,41 @@ static void fifoFree(Fifo *f) { ...@@ -960,19 +962,41 @@ static void fifoFree(Fifo *f) {
free(f); free(f);
} }
static void _particleProblemSolveContacts(ParticleProblem *p, double dt, double tol) { static int _particleProblemSolveContacts(ParticleProblem *p, double dt, double tol) {
int converged = 0; int converged = 0;
int iter = 0;
while (! converged) { 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; converged = 1;
for (int ic = vectorSize(p->contacts)-1; ic >= 0; --ic) { for (int ic = vectorSize(p->contacts)-1; ic >= 0; --ic) {
Contact *c = &p->contacts[ic]; Contact *c = &p->contacts[ic];
//if (iter%2 == 0)
//c = &p->contacts[vectorSize(p->contacts)-1-ic];
int conv; int conv;
switch (c->type) { switch (c->type) {
case PARTICLE_PARTICLE : case PARTICLE_PARTICLE :
conv = contactParticleParticleSolve(c, p, dt,tol); conv = contactParticleParticleSolve(c, p, dt,tol);
break; break;
case PARTICLE_DISK : case PARTICLE_DISK :
conv = contactParticleDiskSolve(c, p, dt, tol); //conv = contactParticleDiskSolve(c, p, dt, tol);
break; break;
case PARTICLE_SEGMENT : case PARTICLE_SEGMENT :
conv = contactParticleSegmentSolve(c, p, dt, tol); conv = contactParticleSegmentSolve(c, p, dt, tol);
...@@ -988,6 +1012,7 @@ static void _particleProblemSolveContacts(ParticleProblem *p, double dt, double ...@@ -988,6 +1012,7 @@ static void _particleProblemSolveContacts(ParticleProblem *p, double dt, double
} }
} }
} }
return 0;
} }
static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, double tol) { static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, double tol) {
...@@ -1079,21 +1104,41 @@ double particleProblemMaxDt(const ParticleProblem *p, double alert) { ...@@ -1079,21 +1104,41 @@ double particleProblemMaxDt(const ParticleProblem *p, double alert) {
return alert / q; 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); _particleProblemGenContacts(p, alert);
if (p->use_queue) /*if (p->use_queue)
_particleProblemSolveContactsQueue(p,dt,tol); _particleProblemSolveContactsQueue(p,dt,tol);
else else*/
_particleProblemSolveContacts(p,dt,tol); return _particleProblemSolveContacts(p,dt,tol);
} }
void particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit) static void particleProblemApplyCtTranslation(ParticleProblem *p)
{
#if FRICTION_ENABLED
for (int i = 0; i < vectorSize(p->contacts); ++i) {
const Contact *c = p->contacts+i;
double t[2] = {c->n[1], -c->n[0]};
double ct = c->ct;
if(c->type == PARTICLE_PARTICLE) {
coordAdd(&p->velocity[c->o0 * DIMENSION], -ct*c->a0, t);
coordAdd(&p->velocity[c->o0 * DIMENSION], ct*c->a1, t);
}
else {
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*c->a1, t);
}
}
#endif
}
int 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 j = 0; j < vectorSize(p->particles); ++j)
for (size_t i = 0; i < DIMENSION; ++i) for (size_t i = 0; i < DIMENSION; ++i)
p->velocity[j * DIMENSION + i] += p->externalForces[j * DIMENSION + i] * dt / p->particles[j].m; p->velocity[j * DIMENSION + i] += p->externalForces[j * DIMENSION + i] * dt / p->particles[j].m;
particleProblemSolve(p, alert, dt, tol, maxit); r = particleProblemSolve(p, alert, dt, tol, maxit);
for (size_t i = 0; i < vectorSize(p->position); ++i){ for (size_t i = 0; i < vectorSize(p->position); ++i){
p->position[i] += p->velocity[i] * dt; p->position[i] += p->velocity[i] * dt;
} }
...@@ -1114,6 +1159,9 @@ void particleProblemIterate(ParticleProblem *p, double alert, double dt, double ...@@ -1114,6 +1159,9 @@ void particleProblemIterate(ParticleProblem *p, double alert, double dt, double
} }
} }
#endif #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) { size_t particleProblemAddBoundarySegmentTagId(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], int tag, int materialTag) {
......
...@@ -31,8 +31,8 @@ ParticleProblem *particleProblemNew(); ...@@ -31,8 +31,8 @@ ParticleProblem *particleProblemNew();
void particleProblemDelete(ParticleProblem *p); void particleProblemDelete(ParticleProblem *p);
void particleProblemLoad(ParticleProblem *p, const char *filename); void particleProblemLoad(ParticleProblem *p, const char *filename);
void particleProblemWrite(const 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); int particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit);
void particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit); int particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit);
double particleProblemMaxDt(const ParticleProblem *p, double alert); double particleProblemMaxDt(const ParticleProblem *p, double alert);
void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r, double m, const char *material); 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); size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x0[DIMENSION], double r, const char *tag, const char *material);
......
...@@ -53,7 +53,8 @@ def genInitialPosition(filename, r, rout, rhop) : ...@@ -53,7 +53,8 @@ def genInitialPosition(filename, r, rout, rhop) :
R2 = x**2 + y**2 R2 = x**2 + y**2
for i in range(x.shape[0]) : for i in range(x.shape[0]) :
if y[i]<0.005: if y[i]<0.005:
r1 = np.random.normal(r,0.05*r) #r1 = np.random.normal(r,0.05*r)
r1 = r
#Addition of an particle object at each point #Addition of an particle object at each point
p.add_particle((x[i], y[i]), r1, r1**2 * np.pi * rhop,("sand" if i%2 == 0 else "rock")); p.add_particle((x[i], y[i]), r1, r1**2 * np.pi * rhop,("sand" if i%2 == 0 else "rock"));
p.write_vtk(filename,0,0) p.write_vtk(filename,0,0)
...@@ -63,7 +64,7 @@ ii = 0 ...@@ -63,7 +64,7 @@ ii = 0
#physical parameters #physical parameters
g = -9.81 # gravity g = -9.81 # gravity
rhop = 2800 # grains density rhop = 2800 # grains density
tEnd = 10 # final time tEnd = 1000 # final time
#numerical parameters #numerical parameters
dt = 1e-3 # time step dt = 1e-3 # time step
...@@ -86,14 +87,14 @@ else : ...@@ -86,14 +87,14 @@ else :
p = scontact.ParticleProblem(2,True) p = scontact.ParticleProblem(2,True)
p.read_vtk(outputdir,0) p.read_vtk(outputdir,0)
outf = 5 # number of iterations between output files outf = 1 # number of iterations between output files
#Taking friction into account #Taking friction into account
p.set_friction_coefficient(1.5,"sand","rock") p.set_friction_coefficient(1.5,"sand","rock")
p.set_friction_coefficient(1.5,"sand","sand") p.set_friction_coefficient(1.5,"sand","sand")
p.set_friction_coefficient(1.5,"rock","rock") p.set_friction_coefficient(1.5,"rock","rock")
p.set_friction_relaxation(1.0) p.set_friction_relaxation(1.)
p.set_tangent_boundary_velocity("Outer",v) p.set_tangent_boundary_velocity("Outer",v)
print(p.velocity().shape[0]) print(p.segments()[:,:4])
#Computation loop #Computation loop
forces = np.zeros_like(p.velocity()) forces = np.zeros_like(p.velocity())
k=0 k=0
...@@ -111,7 +112,11 @@ while t < tEnd : ...@@ -111,7 +112,11 @@ while t < tEnd :
k = k+1 k = k+1
for i in range(nsub) : for i in range(nsub) :
tol = 1e-6 tol = 1e-6
p.iterate(dt/nsub, forces, tol) r = p.iterate(dt/nsub, forces, tol)
if r==1 :
ioutput = int(ii/outf) + 2
p.write_vtk(outputdir, ioutput, t)
exit(0)
t += dt t += dt
#Output files writting #Output files writting
if ii %outf == 0 : if ii %outf == 0 :
......
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