Commit d4527ca1 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

preparing for release

parents aa16cc8f 0c5720b6
......@@ -1728,7 +1728,7 @@ void fluid_problem_after_import(FluidProblem *problem)
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw);
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact, long int *elid, int init, int reload)
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *state, double *contact, long int *elid, int init, int reload)
{
struct timespec tic, toc;
......@@ -1761,16 +1761,17 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->particle_element_id[i]=elid[i];
}
}
mesh_tree_particle_to_mesh(problem->mesh_tree, n, position, problem->particle_element_id, problem->particle_uvw);
int rowl = (D==2 ? 5 : 9);
for (int i = 0; i < n; ++i) {
problem->particle_mass[i] = mass[i];
problem->particle_volume[i] = volume[i];
for (int k = 0; k < D; ++k) {
problem->particle_position[i*D+k] = position[i*D+k];
problem->particle_velocity[i*D+k] = velocity[i*D+k];
problem->particle_position[i*D+k] = state[i*rowl+k];
problem->particle_velocity[i*D+k] = state[i*rowl+D+k];
problem->particle_contact[i*D+k] = contact[i*D+k];
}
}
mesh_tree_particle_to_mesh(problem->mesh_tree, n, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
if (!reload){
for (int i=0; i<problem->mesh->n_nodes; ++i){
problem->oldporosity[i] = problem->porosity[i];
......
......@@ -104,7 +104,7 @@ struct FluidProblem {
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact, long int *elid, int init, int reload);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *state, double *contact, long int *elid, int init, int reload);
void fluid_problem_set_bulk_force(FluidProblem *problem, double *force);
void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
void fluid_problem_free(FluidProblem *problem);
......
......@@ -390,8 +390,8 @@ class FluidProblem :
for i in range(nsub) :
self._lib.fluid_problem_advance_concentration(self._fp,c_double(dt/nsub))
def set_particles(self, mass, volume, position, velocity, contact, init=False, reload = False):
"""Sets the location of the grains in the mesh and compute the porosity in each cell.
def set_particles(self, mass, volume, state, contact, init=False, reload = False):
"""Set location of the grains in the mesh and compute the porosity in each cell.
Keyword arguments:
mass -- List of particles mass
......@@ -402,7 +402,7 @@ class FluidProblem :
reload -- Optional arguments specifying if the simulation is reloaded or if it is the first iteration
"""
self.n_particles = mass.shape[0]
self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),_np2c(mass),_np2c(volume),_np2c(position),_np2c(velocity),_np2c(contact),None,c_int(init),c_int(reload))
self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),_np2c(mass),_np2c(volume),_np2c(state),_np2c(contact),None,c_int(init),c_int(reload))
def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
......
......@@ -58,6 +58,28 @@ def _np2c(a,dtype=np.float64) :
r.tmp = tmp
return r
class ParticleState(np.ndarray) :
@property
def x(self) :
dim = 2 if self.shape[1] == 5 else 3
return self[:,:dim]
@property
def v(self) :
dim = 2 if self.shape[1] == 5 else 3
return self[:,dim:2*dim]
@property
def omega(self) :
dim = 2 if self.shape[1] == 5 else 3
return self[:,2*dim:]
@classmethod
def new(cls, nparticles,dim) :
state = cls((nparticles, 5 if dim == 2 else 9))
return state
class ParticleProblem :
"""Creates the numerical structure containing all the physical particles that appear in the problem"""
......@@ -152,18 +174,21 @@ class ParticleProblem :
"""Returns the list of particle masses."""
return self._get_matrix("Particle", 2)[:, 1][:,None]
def position(self):
"""Returns the list of particle centre position vectors."""
return self._get_matrix("Position", self._dim)
def velocity(self):
"""Returns the list of particle velocity vectors."""
return self._get_matrix("Velocity", self._dim)
def state(self) :
state = ParticleState.new(self.n_particles(),self._dim)
state.x[:] = self._get_matrix("Position", self._dim)
state.v[:] = self._get_matrix("Velocity", self._dim)
if self._friction_enabled and self._rotation_enabled :
state.omega[:] = self._get_matrix("Omega", 1 if self._dim==2 else 3)
else :
state.omega[:] = 0
return state
def omega(self):
"""Returns the list of particle angular velocity."""
assert self._friction_enabled and self._rotation_enabled
return (self._get_matrix("Omega", 1)) if self._dim==2 else self._get_matrix("Omega", 3)
def set_state(self, state) :
self._get_matrix("Position", self._dim)[:] = state.x
self._get_matrix("Velocity", self._dim)[:] = state.v
if self._friction_enabled and self._rotation_enabled :
self._get_matrix("Omega", 1 if self._dim==2 else 3)[:] = state.omega
def particle_material(self):
"""Returns the list of particle materials."""
......@@ -265,7 +290,6 @@ class ParticleProblem :
disks["vt"][disks["tag"]==tag] = vt
segments["vt"][segments["tag"]==tag] = vt
def set_friction_coefficient(self, mu, mat0="default",mat1="default") :
""" Sets the friction coefficient between two materials.
......@@ -321,13 +345,10 @@ class ParticleProblem :
i -- Number of the fiel to write
t -- Time at which the simulation is
"""
v = self.velocity()
if self._rotation_enabled:
omega = self.omega()
else :
omega = np.zeros((v.shape[0],1 if self._dim==2 else 3))
x = self.position()
state = self.state()
v = state.v
omega = state.omega
x = state.x
material = self._get_idx("ParticleMaterial").reshape(-1,1)
forced = self._get_idx("ForcedFlag").reshape(-1,1)
if self._dim == 2 :
......@@ -582,3 +603,6 @@ class ParticleProblem :
(v1[0]+shift[0],v1[1]+shift[1],v1[2]+shift[2]),
(v2[0]+shift[0],v2[1]+shift[1],v2[2]+shift[2]),
stag,material)
def dim(self) :
return self._dim
......@@ -34,11 +34,12 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_
max_split -- maximum number of times the time step can be further split if conergence is not reached
"""
# Compute free velocities of the grains
state = particles.state()
if f is None:
f = np.zeros_like(particles.velocity())
f = np.zeros_like(state.v)
# Estimation of the solid sub-time steps
if particles.r() is not None and particles.r().size != 0:
vn = particles.velocity() + f*dt / particles.mass()
vn = state.v + f*dt / particles.mass()
vmax = np.max(np.linalg.norm(vn,axis=1))
nsub = max(min_nsub, int(np.ceil((vmax * dt * 8)/min(particles.r()))))
else:
......@@ -77,13 +78,12 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
after_sub_iter -- callback to execute once a sub iteration has been made
max_split -- maximum number of times the time step can be further split if conergence is not reached
"""
state0 = particles.state()
if external_particles_forces is None:
external_particles_forces = np.zeros_like(particles.velocity())
external_particles_forces = np.zeros_like(state0.v)
# Copy the fluid solution of the previous time step before computing the prediction
sf0 = np.copy(fluid.solution())
# Copy the solid solution of the previous time step before computing the prediction
vp0 = np.copy(particles.velocity())
pp0 = np.copy(particles.position())
cp0 = np.copy(particles.contact_forces())
if particles._rotation_enabled :
op0 = np.copy(particles.omega())
......@@ -95,10 +95,12 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
sf1 = np.copy(fluid.solution())
f0 = np.copy(fluid.compute_node_force(dt))+external_particles_forces
_advance_particles(particles,f0,dt,min_nsub,contact_tol,max_split=max_split)
state1 = particles.state()
state1.x[:] = state0.x
particles.set_state(state1)
particles.contact_forces()[:] = cp0
# Keep same contact forces and positions while giving to the fluid the predicted solid velocities to compute the correction
particles.position()[:,:] = pp0
particles.contact_forces()[:,:] = cp0
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(), particles.contact_forces(), reload=True)
fluid.set_particles(particles.mass(), particles.volume(), particles.state(), particles.contact_forces(), reload=True)
# Corrector
#
......@@ -114,12 +116,10 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
if fluid.n_fluids() == 2 :
fluid.advance_concentration(dt)
# Reset solid velocities
particles.velocity()[:,:] = vp0
if particles._rotation_enabled :
particles.omega()[:,:] = op0
particles.set_state(state0)
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_split=max_split)
# Give to the fluid the solid information
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(), particles.contact_forces())
fluid.set_particles(particles.mass(), particles.volume(), particles.state(), particles.contact_forces())
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None,max_split=1) :
"""Migflow solver: the solver type depends on the fluid and particles given as arguments.
......@@ -142,7 +142,7 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
if particles is None and fluid is None:
raise ValueError("fluid and particles are both None: not possible to iterate!")
if particles is not None and external_particles_forces is not None:
if external_particles_forces.shape!=particles.velocity().shape:
if external_particles_forces.shape!=(particles.n_particles(),particles.dim()):
raise ValueError("external_particles_forces must have shape (number of particles,dimension!")
if (particles is None or particles.r() is None or particles.r().size == 0) and fluid is not None:
fluid.implicit_euler(dt)
......@@ -155,14 +155,14 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
if fluid is not None and particles is not None and particles.r() is not None and particles.r().size !=0:
if fixed_grains:
fluid.implicit_euler(dt)
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(),-fluid.compute_node_force(dt))
fluid.set_particles(particles.mass(), particles.volume(), particles.state(),-fluid.compute_node_force(dt))
else:
predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter,max_split=max_split)
def particle_changed(fluid,p) :
co = fluid.old_porosity().copy()
c = fluid.porosity().copy()
fluid.set_particles(p.mass(),p.volume(), p.position(), p.velocity(), p.contact_forces())
fluid.set_particles(p.mass(),p.volume(), p.state(), p.contact_forces())
c2 = fluid.porosity()
fluid.old_porosity()[:] = co+c2-c
fluid.velocity()[:] *= c2/c
......
......@@ -82,13 +82,13 @@ fluid.set_wall_boundary("Left")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Right")
# Set locations of the grains in the mesh and compute the porosity in each computation cell.
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces(),init=True)
# Write output files for post-visualization.
fluid.export_vtk(outputdir,0,0)
tic = time.time()
G = np.zeros_like(p.velocity())
G = np.zeros((p.n_particles(),p.dim())
G[:,1] = g[1]*p.mass()[:,0]
#
# COMPUTATION LOOP
......
......@@ -113,7 +113,7 @@ fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Left",0,0)
fluid.set_strong_boundary("Right",0,0)
# Set locations of the grains in the mesh and compute the porosity in each computation cell.
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces(), init=True)
fluid.set_particles(p.mass(), p.volume(), p.state(), p.contact_forces(), init=True)
# Pressure initialization
fluid.solution()[:,2] = (H-fluid.coordinates()[:,1])*(-g[1])*rhof
# Write output files for post-visualization.
......@@ -121,7 +121,7 @@ fluid.export_vtk(outputdir,0,0)
tic = time.time()
G = np.zeros_like(p.velocity())
G = np.zeros((p.n_particles(),p.dim())
G[:,1] = g[1]*p.mass()[:,0]
#
# COMPUTATION LOOP
......
......@@ -130,7 +130,7 @@ fluid.set_wall_boundary("Left")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Right")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces(), init=True)
fluid.set_particles(p.mass(), p.volume(), p.state(), p.contact_forces(), init=True)
#Access to fluid fields and node coordiantes
s = fluid.solution()
x = fluid.coordinates()
......@@ -142,7 +142,7 @@ fluid.set_concentration_cg(c)
fluid.export_vtk(outputdir,0,0)
tic = time.time()
G = np.zeros_like(p.velocity())
G = np.zeros((p.n_particles(),p.dim())
G[:,1] = g[1]*p.mass()[:,0]
#
# COMPUTATION LOOP
......
......@@ -83,7 +83,7 @@ fluid.set_wall_boundary("Left")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Right")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces(),init=True)
# Access to fluid fields and node coordiantes
s = fluid.solution()
x = fluid.coordinates()
......@@ -96,7 +96,7 @@ fluid.set_concentration_cg(c)
fluid.export_vtk(outputdir,0,0)
tic = time.time()
G = np.zeros_like(p.velocity())
G = np.zeros((p.n_particles(),p.dim())
G[:,1] = g[1]*p.mass()[:,0]
#
# COMPUTATION LOOP
......
......@@ -68,7 +68,7 @@ p.set_friction_coefficient(0.9,"Sand","Sand")#Particle-Particle
t = 0
ii = 0
outf = 5
G = np.zeros_like(p.velocity())
G = np.zeros((p.n_particles(),p.dim())
G[:,1] = g[1]*p.mass()[:,0]
#
# COMPUTATION LOOP
......
......@@ -93,7 +93,7 @@ print("number of grains = %d"%len(p.r()))
outf = 2000 #number of iterations between output files
tic = time.time()
G = np.zeros_like(p.velocity())
G = np.zeros((p.n_particles(),p.dim())
G[:,1] = p.mass()[:,0]*g-p.velocity()[:,1]
#
# COMPUTATION LOOP
......
......@@ -117,7 +117,7 @@ fluid.set_wall_boundary("Outer",velocity=[0,0],pressure=0)
fluid.set_strong_boundary("Inner",0,lambda x : (x[:,1]/rout)*0.1)
fluid.set_strong_boundary("Inner",1,lambda x : (-x[:,0]/rout)*0.1)
# Set location of the grains in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces(), init=True)
fluid.set_particles(p.mass(), p.volume(), p.state(), p.contact_forces(), init=True)
fluid.export_vtk(outputdir,0,0)
tic = time.time()
#
......
......@@ -109,13 +109,13 @@ fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
fluid.set_wall_boundary("Top",pressure=0)
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.set_particles(p.mass(), p.volume(), p.state(), p.contact_forces(),init=True)
fluid.export_vtk(outputdir,0,0)
p.write_vtk(outputdir, 0, 0)
tic = time.time()
G = np.zeros_like(p.velocity())
G = np.zeros((p.n_particles(),2))
G[:,1] = p.mass()[:,0]*g[1]
#
# COMPUTATION LOOP
......
......@@ -104,12 +104,12 @@ fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("X")
fluid.set_wall_boundary("Z")
#Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces(),init=True)
tic = time.time()
fluid.export_vtk(outputdir,0,0)
G = np.zeros_like(p.velocity())
G = np.zeros((p.n_particles(),p.dim())
G[:,1] = p.mass()[:,0]*g[1]
#
# COMPUTATION LOOP
......
......@@ -80,7 +80,7 @@ ii = 0
mv.velocity = -1/tEnd
p.setBoundaryVelocity("Pile", [0, mv.velocity])
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.state())
p.write_vtk(outputdir, ii, t)
p.write_boundary_vtk(outputdir, ii, t)
fluid.write_solution(outputdir, ii, t, "vtk")
......@@ -91,7 +91,7 @@ while t < tEnd :
forces = fluid.solve()
p.iterate(dt, forces)
fluid.set_mesh_position(fluid.mesh_position() + mv.call(fluid.mesh_position()) * dt)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.state())
mv.position += dt * mv.velocity
t += dt
ii += 1
......
......@@ -83,7 +83,7 @@ mv.position = 0.5
t = 0
ii = 0
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.state())
p.write(outputdir, ii, t)
fluid.write_solution(outputdir, ii, t)
fluid.write_position(outputdir, ii, t)
......@@ -96,7 +96,7 @@ while t < tEnd :
forces = fluid.solve()
p.iterate(dt, forces)
fluid.set_mesh_position(fluid.mesh_position() + mv.cb(fluid.mesh_position()) * dt)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.state())
mv.position += dt * mv.velocity
t += dt
ii += 1
......
......@@ -123,7 +123,7 @@ fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces(),init=True)
fluid.solution()[:,2] = (H-(fluid.coordinates()[:,1]+0.25))*(-g[1])*rho
t = 0
......
......@@ -123,7 +123,7 @@ fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces(),init=True)
fluid.solution()[:,2] = (H-(fluid.coordinates()[:,1]+0.25))*(-g[1])*rho
t = 0
......
......@@ -116,7 +116,7 @@ fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces(),init=True)
t = 0
fluid.solution()[:,2] = (fluid.coordinates()[:,1]-0.6)*rho*g[1]
......
......@@ -82,7 +82,7 @@ fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
fluid.import_vtk(outputrel+"/fluid_%05d.vtu"%iReload)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),reload=1)
fluid.set_particles(p.mass(), p.volume(), p.state(),reload=1)
ii = (iReload-1)*outf+1
t=ii*dt
......@@ -111,7 +111,7 @@ while t < tEnd :
p.iterate(dt/nsub, forces)
t += dt
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.state())
# Output files writting
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