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

wip

parent 2fef2d8b
Pipeline #8677 failed with stages
in 2 minutes and 1 second
......@@ -1234,7 +1234,7 @@ static void fluid_problem_compute_node_volume(FluidProblem *problem) {
free(periodic_node_volume);
}
static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity, int n_particles, double *particle_position, double *particle_volume, int *particle_element_id, double *particle_uvw)
static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity, int n_particles, double *particle_position, double *particle_volume, int *particle_element_id, double *particle_uvw, int *flag)
{
int independent_nodes = mesh->n_nodes - mesh->n_periodic;
double* periodic_porosity = malloc(sizeof(double)*independent_nodes);
......@@ -1244,7 +1244,7 @@ static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity,
int* periodic_mapping = mesh->periodic_mapping;
for (int i = 0; i < n_particles; ++i) {
if(particle_element_id[i] == -1) continue;
if(particle_element_id[i] == -1 || (flag && !flag[i])) continue;
double sf[N_SF];
shape_functions(&particle_uvw[i*D],sf);
const int *el = &mesh->elements[particle_element_id[i]*N_N];
......@@ -1715,10 +1715,10 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
}
double *old_particle_uvw = malloc(sizeof(double)*D*old_n_particles);
mesh_tree_particle_to_mesh(problem->mesh_tree, old_n_particles, old_particle_position, old_particle_element_id, old_particle_uvw);
compute_porosity(problem->mesh, problem->node_volume, problem->oldporosity, old_n_particles, old_particle_position, old_particle_volume, old_particle_element_id, old_particle_uvw);
compute_porosity(problem->mesh, problem->node_volume, problem->oldporosity, old_n_particles, old_particle_position, old_particle_volume, old_particle_element_id, old_particle_uvw, NULL);
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
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);
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, NULL);
free(old_particle_element_id);
free(old_particle_uvw);
......@@ -1738,21 +1738,15 @@ void fluid_problem_after_import(FluidProblem *problem)
if(problem->mesh_tree)
mesh_tree_free(problem->mesh_tree);
problem->mesh_tree = mesh_tree_create(problem->mesh);
int n_fields = fluid_problem_n_fields(problem);
for (int i = 0; i < problem->n_particles; ++i)
problem->particle_element_id[i] = -1;
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_node_volume(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);
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, NULL);
}
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;
int n_fields = fluid_problem_n_fields(problem);
if(problem->n_particles != n) {
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact) {
if (problem->n_particles != n) {
problem->n_particles = n;
free(problem->particle_uvw);
free(problem->particle_element_id);
......@@ -1764,9 +1758,6 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
free(problem->particle_contact);
problem->particle_uvw = (double*)malloc(sizeof(double)*D*n);
problem->particle_element_id = (int*)malloc(sizeof(int)*n);
for (int i = 0; i < n; ++i) {
problem->particle_element_id[i]=-1;
}
problem->particle_position = (double*)malloc(sizeof(double)*D*n);
problem->particle_mass = (double*)malloc(sizeof(double)*n);
problem->particle_force = (double*)malloc(sizeof(double)*n*D);
......@@ -1774,40 +1765,76 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->particle_velocity = (double*)malloc(sizeof(double)*n*D);
problem->particle_contact = (double*)malloc(sizeof(double)*n*D);
}
if (elid) {
for (int i = 0; i < n; ++i) {
problem->particle_element_id[i]=elid[i];
for (int i = 0; i < n; ++i) {
problem->particle_element_id[i]=-1;
}
mesh_tree_particle_to_mesh(problem->mesh_tree, n, position, problem->particle_element_id, problem->particle_uvw);
double *c = malloc(sizeof(double)*problem->mesh->n_nodes);
for (int i = 0; i < problem->mesh->n_nodes; ++i) {
c[i] = problem->porosity[i];
}
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->n_particles, position, volume, problem->particle_element_id, problem->particle_uvw, NULL);
int n_fields = fluid_problem_n_fields(problem);
for (int i = 0; i < problem->mesh->n_nodes; ++i) {
problem->oldporosity[i] += problem->porosity[i]-c[i];
for (int d = 0; d < DIMENSION; ++d) {
problem->solution[n_fields*i+U+d] *= problem->porosity[i]/c[i];
}
}
free(c);
}
void fluid_problem_move_particles(FluidProblem *problem, int n, double *position, double *velocity, double *contact)
{
if (problem->n_particles != n) {
printf("the number of particle has changed and the \"particle_changed\" function has not been called.\n");
exit(1);
}
for (int i = 0; i < problem->mesh->n_nodes; i++) {
problem->oldporosity[i] = problem->porosity[i];
}
int *old_element_id = malloc(sizeof(int)*problem->n_particles);
for (int i = 0; i < n; ++i) {
old_element_id[i] = problem->particle_element_id[i];
}
mesh_tree_particle_to_mesh(problem->mesh_tree, n, position, problem->particle_element_id, problem->particle_uvw);
int *particle_enter = malloc(sizeof(int)*n);
int *particle_leave = malloc(sizeof(int)*n);
for (int i = 0; i < n; ++i) {
particle_enter[i] = problem->particle_element_id[i] >= 0 && old_element_id[i] < 0;
particle_leave[i] = problem->particle_element_id[i] < 0 && old_element_id[i] >= 0;
}
double *particle_enter_porosity = malloc(sizeof(double)*problem->mesh->n_nodes);
double *particle_leave_porosity = malloc(sizeof(double)*problem->mesh->n_nodes);
compute_porosity(problem->mesh, problem->node_volume, particle_enter_porosity, problem->n_particles, position,
problem->particle_volume, problem->particle_element_id, problem->particle_uvw, particle_enter);
compute_porosity(problem->mesh, problem->node_volume, particle_leave_porosity, problem->n_particles, problem->particle_position,
problem->particle_volume, problem->particle_element_id, problem->particle_uvw, particle_leave);
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->n_particles, position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw, NULL);
for (int i = 0; i < problem->mesh->n_nodes; i++) {
problem->oldporosity[i] += particle_enter_porosity[i]-particle_leave_porosity[i];
}
int nf = fluid_problem_n_fields(problem);
for (int i = 0; i < problem->mesh->n_nodes; i++) {
for (int d = 0; d < DIMENSION; ++d) {
problem->solution[i*nf+U+d] *= problem->porosity[i]/(problem->porosity[i]-particle_enter_porosity[i]+particle_leave_porosity[i]);
}
}
free(particle_enter_porosity);
free(particle_leave_porosity);
free(old_element_id);
free(particle_enter);
free(particle_leave);
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_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];
}
}
if (elid) {
for (int i = 0; i < n; ++i) {
if(elid[i] == -1){
problem->particle_element_id[i]=elid[i];
}
}
}
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);
if (init){
for (int i=0; i<problem->mesh->n_nodes; ++i){
problem->oldporosity[i] = problem->porosity[i];
}
}
}
double *fluid_problem_particle_uvw(FluidProblem *problem) {
......@@ -1824,6 +1851,11 @@ double *fluid_problem_particle_position(FluidProblem *problem)
return problem->particle_position;
}
double *fluid_problem_particle_velocity(FluidProblem *problem)
{
return problem->particle_velocity;
}
double *fluid_problem_particle_volume(FluidProblem *problem)
{
return problem->particle_volume;
......
......@@ -104,7 +104,8 @@ 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 *position, double *velocity, double *contact);
void fluid_problem_move_particles(FluidProblem *problem, int n, double *position, double *velocity, double *contact);
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);
......
......@@ -217,7 +217,7 @@ class FluidProblem :
self.strong_cb_ref.append(bndcb)
self._lib.fluid_problem_set_strong_boundary(self._fp, tag.encode(), c_int(field), c_int(row), bndcb)
def adapt_mesh(self, lcmax, lcmin, n_el, old_n_particles, old_particle_position, old_particle_volume, dt, cmax=1, cmin=1) :
def adapt_mesh(self, lcmax, lcmin, n_el, old_n_particles, old_particle_position, old_particle_volume, cmax=1, cmin=1) :
"""Adapts the mesh.
Keyword arguments:
......@@ -227,11 +227,10 @@ class FluidProblem :
old_n_particles -- Number of particles at the previous time step
old_particle_position -- Position of the particles at the previous time step
old_particle_volume -- Volume of the particles at the previous time step
dt -- time step
cmax -- Optional argument to weigh maximum gradient used in the adaptation formula
cmin -- Optional argument to weigh minimum gradient used in the adaptation formula
"""
self._lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin), c_int(old_n_particles), _np2c(old_particle_position), _np2c(old_particle_volume), c_double(dt))
self._lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin), c_int(old_n_particles), _np2c(old_particle_position), _np2c(old_particle_volume))
self.sys = None
def _mesh_boundaries(self) :
......@@ -278,7 +277,7 @@ class FluidProblem :
offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0])).astype(np.int32)
VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data,cell_data)
def import_vtk(self, filename, dt = 1) :
def import_vtk(self, filename) :
"""Reads output file to reload computation.
Keyword arguments:
......@@ -304,7 +303,7 @@ class FluidProblem :
self.concentration_dg()[:] = cdata["concentration"]
self.porosity()[:] = data["porosity"]
self.old_porosity()[:] = data["old_porosity"]
self._lib.fluid_problem_after_import(self._fp,c_double(dt))
self._lib.fluid_problem_after_import(self._fp)
self.sys = None
def compute_node_force(self, dt) :
......@@ -392,7 +391,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, position, velocity, contact, elid = None, init=False, reload = False):
def set_particles(self, mass, volume, position, velocity, contact):
"""Set location of the grains in the mesh and compute the porosity in each cell.
Keyword arguments:
......@@ -400,16 +399,20 @@ class FluidProblem :
volume -- List of particles volume
position -- List of particles centre positions
velocity -- List of particles velocity
init -- If True sets the old_porosity to zero (has to be used at the first call)
reload -- Optional arguments specifying if the simulation is reloaded or if it is the first iteration
contact -- List of particles contact resultant forces
"""
self.n_particles = mass.shape[0]
elemid = None
if elid is not None:
elemid = _np2c(elid,dtype=int)
self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),_np2c(mass),_np2c(volume),_np2c(position),_np2c(velocity),_np2c(contact),elemid,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))
def move_particles(self, position, velocity, contact):
"""Set location of the grains in the mesh and compute the porosity in each cell.
Keyword arguments:
position -- List of particles centre positions
velocity -- List of particles velocity
contact -- List of particles contact resultant forces
"""
self._lib.fluid_problem_move_particles(self._fp,c_int(velocity.shape[0]),_np2c(position),_np2c(velocity),_np2c(contact))
def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
f = getattr(self._lib,"fluid_problem_"+f_name)
......@@ -505,9 +508,13 @@ class FluidProblem :
return self._get_matrix("particle_uvw",self.n_particles,self._dim)
def particle_position(self) :
"""Gives access to the fluid field value at the mesh nodes."""
"""Gives access to stored paricles position."""
return self._get_matrix("particle_position", self.n_particles, self._dim)
def particle_velocity(self) :
"""Gives access to the stored particle velocity."""
return self._get_matrix("particle_velocity", self.n_particles, self._dim)
def particle_volume(self) :
"""Gives access to the fluid field value at the mesh nodes."""
return self._get_matrix("particle_volume", self.n_particles, 1)
......
......@@ -39,12 +39,12 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_
f = np.zeros_like(v)
# Estimation of the solid sub-time steps
if particles.r() is not None and particles.r().size != 0:
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()))))
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:
nsub = 1
vmax = 0
nsub = 1
vmax = 0
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(particles.r()) if (particles.r() is not None and particles.r().size !=0) else 0 )
nsubt2 = nsubt + [nsub]
#If time step was split too many times, ignore the convergence and move the grains
......@@ -79,7 +79,6 @@ 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
"""
flags_inside = fluid.particle_element_id().copy()
particles.save_state()
if external_particles_forces is None:
external_particles_forces = np.zeros_like(particles.velocity())
......@@ -102,7 +101,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
_advance_particles(particles,f0,dt,min_nsub,contact_tol,max_split=max_split)
# 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(), x0, particles.velocity(), cp0, reload=True)
fluid.particle_velocity()[:] = particles.velocity()
# Corrector
#
......@@ -121,12 +120,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
particles.restore_state()
_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(),elid=flags_inside)
oldporosity1 = fluid.old_porosity().copy()
porosity1 = fluid.porosity().copy()
fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(), particles.contact_forces())
fluid.old_porosity()[:] = oldporosity1 + fluid.porosity() - porosity1
fluid.velocity()[:] *= fluid.porosity()/porosity1
fluid.move_particles(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) :
......@@ -163,17 +157,6 @@ 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.move_particles(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)
def particle_changed(fluid, p, dt) :
""" Modifies the porosity and fluid velocity fields when particles have been removed or added in the fluid domain
"""
c = fluid.porosity().copy()
oldporosity0 = fluid.old_porosity().copy()
fluid.set_particles(p.mass(),p.volume(), p.position(), p.velocity(), p.contact_forces())
c2 = fluid.porosity()
fluid.old_porosity()[:] = oldporosity0 + c2 -c
fluid.velocity()[:] *= c2/c
......@@ -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.position(), p.velocity(), p.contact_forces())
fluid.export_vtk(outputdir,0,0)
tic = time.time()
#
......
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