Skip to content
GitLab
Commit c80fde53 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

correcting bug with writing forces instead of impulses

parent 490f2c2c
Pipeline #8482 failed with stages
in 2 minutes and 10 seconds
Hide whitespace changes
Inline Side-by-side
python/scontact.py
View file @ c80fde53
......@@ -115,7 +115,7 @@ class ParticleProblem :
return np.ctypeslib.as_array(buf)
def __init__(self, dim, friction_enabled=False, rotation_enabled=True) :
"""Builds the particles contanier structure.
"""Builds the particles problem structure.
Keyword arguments:
dim -- Dimension of the problem
......@@ -175,6 +175,7 @@ class ParticleProblem :
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)
......@@ -185,6 +186,7 @@ class ParticleProblem :
return state
def set_state(self, state) :
"""Sets the positions, velocities and angular velocities of the particle problem"""
self._get_matrix("Position", self._dim)[:] = state.x
self._get_matrix("Velocity", self._dim)[:] = state.v
if self._friction_enabled and self._rotation_enabled :
......@@ -381,6 +383,7 @@ class ParticleProblem :
nameb = name.encode()
fdata.append((nameb,v[:,[0]]))
fdata.append((nameb+b"_t",v[:,[1]]))
fdata.append((nameb+b"_dir_n",basis[:,0,:]))
if self._dim == 3:
fdata.append((nameb+b"_s",v[:,[2]]))
fdata.append((nameb+b"_idx",o))
......@@ -431,12 +434,21 @@ class ParticleProblem :
_,_,_,_,fdata = VTK.read(dirname+"/contacts_%05d.vtp"%i)
ks = ("particle_particle","particle_disk","particle_segment","particle_triangle")
v = []
basis = []
v.append(np.vstack([fdata[k] for k in ks]))
v.append(np.vstack([fdata[k+"_t"] for k in ks]))
if self._dim == 3:
basis.append(np.vstack([fdata[k+"_dir_n"] for k in ks]))
if self._dim == 2:
basis.append(np.vstack([fdata[k+"_dir_n"] for k in ks])[:,[1,0]]*np.array([-1, 1]))
else:
dirt = np.vstack([fdata[k+"_dir_n"] for k in ks])[:,[2,0,1]]*np.array([-1, 1, -1]) - np.vstack([fdata[k+"_dir_n"] for k in ks])*np.einsum("ij,ij->i",np.vstack([fdata[k+"_dir_n"] for k in ks]),np.vstack([fdata[k+"_dir_n"] for k in ks])[:,[2,0,1]]*np.array([-1, 1, -1]))[:,np.newaxis]
nnorm = np.linalg.norm(dirt,axis=1)
nnorm[nnorm==0] = 1
basis.append(dirt/nnorm[:,np.newaxis])
basis.append(np.cross(np.vstack([fdata[k+"_dir_n"] for k in ks]),dirt/nnorm[:,np.newaxis]))
v.append(np.vstack([fdata[k+"_s"] for k in ks]))
basis = np.column_stack(basis)
v = np.column_stack(v)*contact_factor
basis = np.zeros((v.shape[0],self._dim,self._dim))
o = np.vstack([fdata[k+"_idx"] for k in ks])
t = np.repeat([0,1,2,3],[fdata[k].shape[0] for k in ks])
self._lib.particleProblemSetContacts(self._p,c_int(t.shape[0]),_np2c(o,np.uint64),_np2c(v),_np2c(basis),_np2c(t,np.int32))
......@@ -481,7 +493,7 @@ class ParticleProblem :
return self._lib.particleProblemAddBoundaryTriangle(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), tag.encode(),material.encode())
def add_particle(self, x, r, m, tag="default",forced = 0):
"""Adds a particle in the particle container.
"""Adds a particle in the particle problem.
Keyword arguments:
x -- Tuple to set the centre particle position
......@@ -493,7 +505,7 @@ class ParticleProblem :
self._lib.particleProblemAddParticle(self._p, self._coord_type(*x), c_double(r), c_double(m), tag.encode(), c_int(forced))
def add_particles(self, x, r, m, v=None, omega=None, tag="default", forced=None, contact_forces=None):
"""Adds particles in the particle container.
"""Adds particles in the particle problem.
Keyword arguments:
x -- Array with the particles's centers position
r -- Array with the particles's radii
......@@ -605,4 +617,5 @@ class ParticleProblem :
stag,material)
def dim(self) :
"""Returns the dimension of the particle problem"""
return self._dim
python/time_integration.py
View file @ c80fde53
......@@ -51,7 +51,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_
for i in range(nsub) :
particles.iterate(dt/nsub, f, tol=contact_tol,force_motion=1)
if after_sub_iter :
after_sub_iter(dt)
after_sub_iter(nsub)
return
# For each sub-time step solve contacts (do NLGS iterations) and move the grains
for i in range(nsub) :
......@@ -61,7 +61,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_
print("Convergence reached for subinvervals of level %d"%(iteration+1))
else :
if after_sub_iter :
after_sub_iter(dt)
after_sub_iter(nsub)
def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5,after_sub_iter=None,max_split=1) :
......@@ -87,6 +87,11 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
cp0 = np.copy(particles.contact_forces())
if particles._rotation_enabled :
op0 = np.copy(particles.omega())
seg0 = np.copy(particles.segments())
disk0 = np.copy(particles.disks())
if particles.dimension == 3 :
tri0 = np.copy(particles.triangles())
# Predictor
#
......@@ -99,6 +104,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
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)
......@@ -117,6 +123,10 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
fluid.advance_concentration(dt)
# Reset solid velocities
particles.set_state(state0)
particles.segments()[:] = seg0
particles.disks()[:] = disk0
if particles.dimension == 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())
......@@ -160,6 +170,8 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
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) :
""" Modifies the porosity and fluid velocity fields when particles have been removed or added in the fluid domain
"""
co = fluid.old_porosity().copy()
c = fluid.porosity().copy()
fluid.set_particles(p.mass(),p.volume(), p.state(), p.contact_forces())
......
scontact/scontact.c
View file @ c80fde53
......@@ -1121,6 +1121,15 @@ int particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol
int particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit,int force_motion, double *externalForces)
{
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
p->contacts[j].dv *= dt;
p->contacts[j].ct *= dt;
#if DIMENSION==3
p->contacts[j].cs *= dt;
#endif
}
double *oldVelocity = malloc(sizeof(double)*vectorSize(p->velocity));
for (size_t j = 0; j < vectorSize(p->particles); ++j){
for (size_t i = 0; i < DIMENSION; ++i){
......@@ -1141,6 +1150,14 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
p->contactForces[j*DIMENSION+i] = 0;
}
}
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
p->contacts[j].dv /= dt;
p->contacts[j].ct /= dt;
#if DIMENSION==3
p->contacts[j].cs /= dt;
#endif
}
free(oldVelocity);
return 0;
}
......@@ -1170,6 +1187,15 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
}
}
free(oldVelocity);
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
p->contacts[j].dv /= dt;
p->contacts[j].ct /= dt;
#if DIMENSION==3
p->contacts[j].cs /= dt;
#endif
}
return 1;
}
......
testcases/drag-2d/drag.py
View file @ c80fde53
......@@ -28,11 +28,11 @@ from migflow import fluid
from migflow import time_integration
import numpy as np
import os
import time
import shutil
import random
import sys
index = int(sys.argv[1])
vitesse = ["36","64","95","130","303","408"]
# Physical parameters
vit = -0.1 #stream velocity
vit = -0.001*int(vitesse[index]) #stream velocity
rhop = 2500 #particles density
rho = 1000 #fluid density
nu = 1e-6 #fluid kinematic viscosity
......@@ -40,11 +40,11 @@ g = np.array([0,-9.81]) #gravity
# Numerical parameters
dt = 1e-3 #time step
t = 0 #initial time
tEnd = 13 #final time
tEnd = 10 #final time
ii = 0 #iteration number
outf = 10 #iterations between data frames
# Define output directory
outputdir = "output"
outputdir = "notOkTrial/" + vitesse[index]
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
filename = outputdir + "/Drag.csv"
......@@ -53,14 +53,13 @@ filename = outputdir + "/Drag.csv"
#
#Injected particles
p = scontact.ParticleProblem(2,True,True)
p.read_vtk("depot/output",10)
state = p.state()
p.remove_particles_flag((state.x[:,1]-p.r()[:,0]>0.5)*(state.x[:,1]+p.r()[:,0] <0.53))
p.read_vtk("depot/outputNotOk",10)
p.remove_particles_flag((p.state().x[:,1]-p.r()[:,0]>0.5)*(p.state().x[:,1]+p.r()[:,0] <0.61))
#Flow particles
p2 = scontact.ParticleProblem(2,True,True)
p2.read_vtk("depot/output",10)
p2.read_vtk("depot/outputNotOk",10)
state2 = p2.state()
p2.remove_particles_flag(state2.x()[:,1]+p2.r()[:,0] <0.5)
p2.remove_particles_flag(state2.x[:,1]+p2.r()[:,0] <0.5)
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner","Left","Right"],material = "Steel")
p2.set_friction_coefficient(0.3,"Sand","Sand")
......@@ -69,38 +68,43 @@ p2.write_vtk(outputdir, 0, 0)
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=0)
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Inner")
fluid.set_wall_boundary("Left", velocity = [0,0])
fluid.set_wall_boundary("Right", velocity = [0,0])
fluid.set_open_boundary("Bottom", velocity = [0,vit])
fluid.set_open_boundary("Top",pressure = 0)
fluid.set_particles(p2.mass(), p2.volume(), state2, p2.contact_forces(),init=True)
fluid.export_vtk(outputdir, 0.,0)
#fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=0)
#fluid.load_msh("mesh.msh")
#fluid.set_wall_boundary("Inner")
#fluid.set_wall_boundary("Left", velocity = [0,0])
#fluid.set_wall_boundary("Right", velocity = [0,0])
#fluid.set_open_boundary("Bottom", velocity = [0,vit])
#fluid.set_open_boundary("Top",pressure = 0)
#fluid.set_particles(p2.mass(), p2.volume(), state2, p2.contact_forces(),init=True)
#fluid.export_vtk(outputdir, 0.,0)
#
# COMPUTATION LOOP
#
def accumulate(F):
F+=np.sum(p2.get_boundary_impulses("Inner"),axis=0)
def accumulate(F,nsub):
F+=np.sum(p2.get_boundary_impulses("Inner"),axis=0)/nsub
while t < tEnd :
state2 = p2.state()
p2.forced_flag()[(state2.x[:,1] + p2.r()[:,0] <= -0.5)] = 1
state2.v[p2.forced_flag()==1,:] = [0, vit]
state2.omega[p2.forced_flag()==1] = 0
if max(state2.x[:,1]+p2.r()[:,0]) < 0.5 and t < 25 :
p2.add_particles(state.x,p.r(),p.mass(),v=state.v,tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.remove_particles_flag( (state2.x[:,1] + p2.r()[:,0] >-0.55))
if t > 2:
state2.v[p2.forced_flag()==1,:] = [0, vit]
else :
state2.v[p2.forced_flag()==1,:] = [0, 0]
state2.omega[p2.forced_flag()==0] = 1
p2.set_state(state2)
if max(p2.state().x[:,1]+p2.r()[:,0]) < 0.5 and t < 9 :
p2.add_particles(p.state().x,p.r(),p.mass(),v=p.state().v,tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.remove_particles_flag( (p2.state().x[:,1] + p2.r()[:,0] >-0.55))
F = np.zeros(2)
time_integration.particle_changed(fluid,p2)
time_integration.iterate(fluid,p2,dt,1,contact_tol=5e-7,external_particles_forces=g*p2.mass(),after_sub_iter=lambda subdt : accumulate(F))
#time_integration.particle_changed(fluid,p2)
time_integration.iterate(None,p2,dt,1,contact_tol=5e-7,external_particles_forces=g*p2.mass(),after_sub_iter=lambda nsub : accumulate(F,nsub))
with open(filename,"a") as file1:
F += fluid.boundary_force()[0,:]
#F += fluid.boundary_force()[0,:]
file1.write( str(F[0])+";"+str(F[1])+";"+str(t)+"\n")
t += dt
# Output files writing
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p2.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir+"/Flow", t, ioutput)
#fluid.export_vtk(outputdir+"/Flow", t, ioutput)
ii += 1
print("%i : %.2g/%.2g" % (ii, t, tEnd))
testcases/drag-3d/drag.py
View file @ c80fde53
......@@ -60,14 +60,14 @@ filename = outputdir + "/Drag.csv"
#Injected particles
p = scontact.ParticleProblem(3,True,True)
p.read_vtk("depot/"+ str(numDep),2)
p.remove_particles_flag((p.state().x[:,1]-p.r()[:,0]>0.5)*(p.state().x[:,1]+p.r()[:,0] <0.52))
state = p.state()
p.remove_particles_flag((state.x[:,1]-p.r()[:,0]>0.5)*(state.x[:,1]+p.r()[:,0] <0.52))
state.v[:,1] = vit
p.set_state(state)
#Flow particles
p2 = scontact.ParticleProblem(3,True,True)
p.read_vtk("depot/"+ str(numDep),2)
state2 = p2.state()
p2.remove_particles_flag(state2.x[:,1]+p2.r()[:,0] <0.5)
p2.remove_particles_flag(p2.state().x[:,1]+p2.r()[:,0] <0.5)
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner","Left", "Right","Front","Rear"],material = "Steel")
p2.set_friction_coefficient(mupart,"Sand","Sand")
......@@ -90,20 +90,21 @@ p2.write_vtk(outputdir, 0, 0)
#
# COMPUTATION LOOP
#
def accumulate(F):
F+=np.sum(p2.get_boundary_impulses("Inner"),axis=0)
def accumulate(F,nsub):
F+=np.sum(p2.get_boundary_impulses("Inner"),axis=0)/nsub
while t < tEnd :
state2 = p2.state()
p2.forced_flag()[(state2.x[:,1] + p2.r()[:,0] <= -0.5)] = 1
state2.v[p2.forced_flag()==1,:] = [0, vit, 0]
state2.omega[p2.forced_flag()==1,:] = [0,0,0]
if max(state2.x[:,1]+p2.r()[:,0]) < 0.5 and t < 21 :
p2.add_particles(state.x,p.r(),p.mass(),v=state.v,tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.remove_particles_flag( (state2.x[:,1] + p2.r()[:,0] >-0.52))
p2.set_state(state2)
if max(p2.state().x[:,1]+p2.r()[:,0]) < 0.5 and t < 21 :
p2.add_particles(p.state().x,p.r(),p.mass(),v=p.state().v,tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.remove_particles_flag( (p2.state().x[:,1] + p2.r()[:,0] >-0.52))
F = np.zeros(3)
#time_integration.particle_changed(fluid,p2)
time_integration.iterate(None,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda subdt : accumulate(F))
time_integration.iterate(None,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda nsub : accumulate(F,nsub))
with open(filename,"a") as file1:
F /= dt
#F += fluid.boundary_force()[0,:]
file1.write( str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+"\n")
t += dt
......
testcases/hopper-2d/hopper.py
View file @ c80fde53
......@@ -81,6 +81,7 @@ p2.remove_particles_flag((state2.x[:,1] + p2.r()[:,0] >1.1))
p2.remove_particles_flag((state2.x[:,0] + p2.r()[:,0] <0.15))
p2.remove_particles_flag((state2.x[:,0] + p2.r()[:,0] >-0.15))
state2.x[:,1] -= 0.6
p2.set_state(state2)
p.clear_boundaries()
p.load_msh_boundaries("mesh.msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh.msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"], material = "Plexi")
......@@ -93,8 +94,8 @@ p.write_vtk(outputdir, 0, 0)
# COMPUTATION LOOP
#
opened = 0
def accumulate(F):
F+=np.sum(p.get_boundary_forces("Inner"),axis=0)
def accumulate(F,nsub):
F+=np.sum(p.get_boundary_forces("Inner"),axis=0)/nsub
while t < tEnd :
if t>=3 and opened == 0:
opened = 1
......@@ -102,13 +103,13 @@ while t < tEnd :
p.load_msh_boundaries("mesh.msh", ["Inner"],material="Steel")
p.load_msh_boundaries("mesh.msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"],material="Plexi")
#Adding new particles
if (max(state.x[:,1])+max(p.r()) <0.5) and t > 3 and t < 25:
p.add_particles(state2.x,p2.r(),p2.mass(),v=state2.v,tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
if (max(p.state().x[:,1])+max(p.r()) <0.5) and t > 3 and t < 25:
p.add_particles(p2.state().x,p2.r(),p2.mass(),v=p2.state().v,tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
#Solving the contacts
F = np.zeros(2)
time_integration.iterate(None,p,dt,min_nsub=1,contact_tol=1e-7,external_particles_forces=g*p.mass(),after_sub_iter= lambda subdt : accumulate(F))
time_integration.iterate(None,p,dt,min_nsub=1,contact_tol=1e-7,external_particles_forces=g*p.mass(),after_sub_iter= lambda nsub : accumulate(F,nsub))
#Removing particles outside the hopper
p.remove_particles_flag( (state.x[:,1] >-0.6))
p.remove_particles_flag( (p.state().x[:,1] >-0.6))
#Computing mean velocity and writing to file
with open(filename,"a") as file1:
file1.write(str(F[0])+";"+str(F[1])+";"
......
testcases/hopper-3d/hopper.py
View file @ c80fde53
......@@ -82,13 +82,13 @@ else :
p2 = scontact.ParticleProblem(3,friction_enabled = True,rotation_enabled=True)
p.read_vtk(outputdir, 0)
p2.read_vtk(outputdir, 0)
state = p.state()
p2.remove_particles_flag((p2.state().x[:,1] + p2.r()[:,0] <1.25))
p2.remove_particles_flag((p2.state().x[:,1] + p2.r()[:,0] >1.1))
p2.remove_particles_flag((p2.state().x[:,0] + p2.r()[:,0] <0.18))
p2.remove_particles_flag((p2.state().x[:,0] + p2.r()[:,0] >-0.18))
state2 = p2.state()
p2.remove_particles_flag((state2.x[:,1] + p2.r()[:,0] <1.25))
p2.remove_particles_flag((state2.x[:,1] + p2.r()[:,0] >1.1))
p2.remove_particles_flag((state2.x[:,0] + p2.r()[:,0] <0.18))
p2.remove_particles_flag((state2.x[:,0] + p2.r()[:,0] >-0.18))
state2.x[:,1] -= 0.5
p2.set_state(state2)
p.clear_boundaries()
p.load_msh_boundaries("mesh" + ".msh", ["Inner","Right","Left","Front","Rear","LockDown"], material = "Steel")
......@@ -106,12 +106,12 @@ while t < tEnd :
p.load_msh_boundaries("mesh" + ".msh", ["Inner","Right","Left","Front","Rear"],material="Steel")
#Adding new particles
if (max(state.x[:,1])+max(p.r()) <0.6) and t > 10 and t < 25:
p.add_particles(state2.x,p2.r(),p2.mass(),v=state2.v,tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
p.add_particles(p2.state().x,p2.r(),p2.mass(),v=p2.state().v,tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
#Solving the contacts
F = np.zeros(3)
time_integration.iterate(None,p,dt,min_nsub=1,contact_tol=1e-6,external_particles_forces=g*p.mass(),after_sub_iter= lambda subdt : accumulate(F))
#Removing particles outside the hopper
p.remove_particles_flag( (state.x[:,1] >-0.6))
p.remove_particles_flag( (p.state().x[:,1] >-0.6))
#Computing mean velocity and writing to file
with open(filename,"a") as file1:
file1.write(str(F[0])+";"+str(F[1])+";"+str(F[2])+";"
......
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