Commit 8668533f authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

replace particles.state() by save_state()/restore_state()

parent cadabba3
Pipeline #8533 failed with stages
in 1 minute and 41 seconds
......@@ -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 *state, double *contact, long int *elid, int init, int reload)
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)
{
struct timespec tic, toc;
......@@ -1761,13 +1761,12 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->particle_element_id[i]=elid[i];
}
}
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] = state[i*rowl+k];
problem->particle_velocity[i*D+k] = state[i*rowl+D+k];
problem->particle_position[i*D+k] = position[i*D+k];
problem->particle_velocity[i*D+k] = velocity[i*D+k];
problem->particle_contact[i*D+k] = contact[i*D+k];
}
}
......
......@@ -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 *state, double *contact, long int *elid, int init, int reload);
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_bulk_force(FluidProblem *problem, double *force);
void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
void fluid_problem_free(FluidProblem *problem);
......
......@@ -390,7 +390,7 @@ 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, state, contact, init=False, reload = False):
def set_particles(self, mass, volume, position, velocity, contact, init=False, reload = False):
"""Set location of the grains in the mesh and compute the porosity in each cell.
Keyword arguments:
......@@ -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(state),_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(position),_np2c(velocity),_np2c(contact),None,c_int(init),c_int(reload))
def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
......
......@@ -58,28 +58,6 @@ 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"""
......@@ -174,16 +152,18 @@ class ParticleProblem :
"""Returns the list of particle masses."""
return self._get_matrix("Particle", 2)[:, 1][:,None]
def state(self) :
"""Returns the positions, velocities and angular velocities of the particle problem"""
state = ParticleState.new(self.n_particles(),self._dim)
state.x[:] = self._get_matrix("Position", self._dim)
state.v[:] = self._get_matrix("Velocity", self._dim)
def position(self):
return self._get_matrix("Position", self._dim)
def velocity(self):
return self._get_matrix("Velocity", self._dim)
def omega(self) :
d = 1 if self._dim==2 else 3
if self._friction_enabled and self._rotation_enabled :
state.omega[:] = self._get_matrix("Omega", 1 if self._dim==2 else 3)
return(self._get_matrix("Omega", d))
else :
state.omega[:] = 0
return state
return np.zeros((self.n_particles(),d))
def set_state(self, state) :
"""Sets the positions, velocities and angular velocities of the particle problem"""
......@@ -192,6 +172,18 @@ class ParticleProblem :
if self._friction_enabled and self._rotation_enabled :
self._get_matrix("Omega", 1 if self._dim==2 else 3)[:] = state.omega
def save_state(self) :
self._saved_velocity = np.copy(self.velocity())
self._saved_position = np.copy(self.position())
if self._friction_enabled and self._rotation_enabled :
self._saved_omega = np.copy(self.omega())
def restore_state(self) :
self.velocity()[:] = self._saved_velocity
self.position()[:] = self._saved_position
if self._friction_enabled and self._rotation_enabled :
self.omega()[:] = self._saved_omega
def particle_material(self):
"""Returns the list of particle materials."""
return self._get_idx("ParticleMaterial")
......
......@@ -34,12 +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()
v = particles.velocity()
if f is None:
f = np.zeros_like(state.v)
f = np.zeros_like(v)
# Estimation of the solid sub-time steps
if particles.r() is not None and particles.r().size != 0:
vn = state.v + f*dt / particles.mass()
vn = 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:
......@@ -78,12 +78,13 @@ 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()
particles.save_state()
if external_particles_forces is None:
external_particles_forces = np.zeros_like(state0.v)
external_particles_forces = np.zeros_like(particles.velocity())
# 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
x0 = np.copy(particles.position())
cp0 = np.copy(particles.contact_forces())
seg0 = np.copy(particles.segments())
disk0 = np.copy(particles.disks())
......@@ -98,13 +99,10 @@ 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
fluid.set_particles(particles.mass(), particles.volume(), particles.state(), particles.contact_forces(), reload=True)
fluid.set_particles(particles.mass(), particles.volume(), x0, particles.velocity(), particles.contact_forces(), reload=True)
# Corrector
#
......@@ -120,14 +118,14 @@ 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.set_state(state0)
particles.restore_state()
particles.segments()[:] = seg0
particles.disks()[:] = disk0
if particles.dim() == 3 :
particles.triangles()[:] = tri0
_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.state(), particles.contact_forces())
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(), 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.
......@@ -163,7 +161,7 @@ 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.state(),-fluid.compute_node_force(dt))
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(),-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)
......@@ -172,7 +170,7 @@ def particle_changed(fluid,p) :
"""
co = fluid.old_porosity().copy()
c = fluid.porosity().copy()
fluid.set_particles(p.mass(),p.volume(), p.state(), p.contact_forces())
fluid.set_particles(p.mass(),p.volume(), p.position(), p.velocity(), p.contact_forces())
c2 = fluid.porosity()
fluid.old_porosity()[:] = co+c2-c
fluid.velocity()[:] *= c2/c
......
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