Skip to content
GitLab
Commit fdce8d90 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

test cases with comments

parent bfb9e45e
Hide whitespace changes
Inline Side-by-side
testcases/couette-2d/betonlmgc.py
View file @ fdce8d90
#!/usr/bin/env python
#variable to use lmgc or not
use_lmgc = 1
......@@ -17,28 +18,38 @@ outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#grains creation and placement + physical boundaries definition
'''r is the radius of the grains
rout is the outer radius of the geometry
rhop is the particles density'''
def genInitialPosition(filename, r, rout, rhop) :
p = scontact2.ParticleProblem()
#Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
#Definition of the points where the grains are located
x = np.arange(rout, -rout, 2 * -r)
y = np.arange(-rout, -2*rout/3., 2 * r)
x, y = np.meshgrid(x, y)
R2 = x**2 + y**2
#condition to be inside the outer boundary
keep = R2 < (rout - r)**2
x = x[keep]
y = y[keep]
#First layer of grains
for i in range(x.shape[0]) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
x = np.arange(rout, -rout, 2 * -r)
y = np.arange(-2*rout/3.+r, -rout/2., 2 * r)
x = np.arange(rout, -rout, -r)
y = np.arange(-2*rout/3.+r, -rout/2., r)
x, y = np.meshgrid(x, y)
R2 = x**2 + y**2
keep = R2 < (rout - r)**2
x = x[keep]
y = y[keep]
#Second layer of grains
for i in range(x.shape[0]) :
p.add_particle((x[i], y[i]), r/2., r**2/4. * np.pi * rhop);
x = np.arange(rout, -rout, 2 * -r/3)
y = np.arange(-rout/2.+r, -5*rout/12., 2 * r/3)
x, y = np.meshgrid(x, y)
......@@ -46,6 +57,7 @@ def genInitialPosition(filename, r, rout, rhop) :
keep = R2 < (rout - r/3)**2
x = x[keep]
y = y[keep]
#Third layer of grains
for i in range(x.shape[0]) :
p.add_particle((x[i], y[i]), r/3, r**2 /9 * np.pi * rhop);
p.write_vtk(filename,0,0)
......@@ -55,61 +67,64 @@ t = 0
ii = 0
#physical parameters
g = -9.81
rho = 1e3
rhop = 1600#1.18e3
nu = 1e-4
tEnd = 100
g = -9.81 # gravity
rho = 1e3 # fluid density
rhop = 1600 # grains density
nu = 1e-4 # kinematic viscosity
tEnd = 100 # final time
#numerical parameters
h = 0.002 # approx r*100 but should match the mesh size
dt = 5e-3#h/V*0.1
epsilon = 1e-4 # ?? not sure ??import pyfefp
h = 0.002 # mesh size
dt = 5e-3 # time step
epsilon = 1e-4 # stabilization parametre
rout = 0.0254
rin = 0.0064
r = 397e-6/2
compacity3d = 0.18
N = int((rout**2 - rin**2) * compacity3d * 3./2 / r**2)
#Geometrical parametres
rout = 0.0254 # outer boundary
r = 397e-6/2 # grains radius
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#Object particles creation
genInitialPosition(outputdir, r, rout, rhop)
if use_lmgc :
friction=0.1
friction=0.1 # friction coefficient
lmgc2Interface.scontactToLmgc90(outputdir,0,friction)
p = lmgc2Interface.LmgcInterface()
else :
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
outf = 5
outf = 5 # number of iterations between output files
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Outer",0,lambda x : -4*x[:, 1]),("Outer",1,lambda x : 4*x[:, 0]),("Inner",0,0.),("Inner",1,0.),("PtFix",2,0.)]
#g is define in fluid_problem
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)
tic = time.clock()
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
#Computation of the new velocities
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
#number of sub time step
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
#Contact solver
for i in range(nsub) :
p.iterate(dt/nsub, forces)
tol = 1e-6
p.iterate(dt/nsub, forces, tol)
#Localisation of the grains in the fluid
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.implicit_euler(dt)
t += dt
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir,ioutput,t)
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
if(ii == 100):
exit()
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
testcases/couette-2d/melangeur.py
View file @ fdce8d90
......@@ -6,26 +6,39 @@ import numpy as np
import os
import time
import shutil
import random
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#grains creation and placement + physical boundaries definition
'''r is the radius of the grains
rout is the outer radius of the geometry
rin is the inner radius of the geometry
rhop is the particles density'''
def genInitialPosition(filename, r, rout, rin, rhop) :
p = scontact2.ParticleProblem()
#Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
#Definition of the points where the grains are located
x = np.arange(rout, -rout, 2 * -r)
x, y = np.meshgrid(x, x)
R2 = x**2 + y**2
#condition to be inside the outer boundary
keep = R2 < (rout - r)**2
x = x[keep]
y = y[keep]
R2 = x**2 + y**2
#condition to be outside the inner boundary
keep = R2 > (rin + r)**2
x = x[keep]
y = y[keep]
for i in range(x.shape[0]) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
if y[i]<0:
rhop1 = random.choice([rhop*.9,1.1*rhop,rhop])
#Addition of an particle object at each point
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1);
p.write_vtk(filename,0,0)
......@@ -33,50 +46,58 @@ t = 0
ii = 0
#physical parameters
g = 0
rho = 1.253e3
rhop = 1000#1.18e3
nu = 1e-4
tEnd = 50
g = 0 # gravity
rho = 1.253e3 # fluid density
rhop = 1000 # grains density
nu = 1e-4 # kinematic viscosity
tEnd = 50 # final time
#numerical parameters
h = 0.002 # approx r*100 but should match the mesh size
dt = 1e-2#h/V*0.1
epsilon = 1e-5 # ?? not sure ??import pyfefp
h = 0.002 # mesh size
dt = 1e-2 # time step
epsilon = 1e-5 # stabilisation parametre
rout = 0.0254
rin = 0.0064
r = 397e-6/2
compacity3d = 0.18
N = int((rout**2 - rin**2) * compacity3d * 3./2 / r**2)
p = scontact2.ParticleProblem()
#geometry parameters
rout = 0.0254 # outer radius
rin = 0.0064 # inner radius
r = 397e-6/2 # grains radius
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#Object particles creation
genInitialPosition(outputdir, r, rout, rin, rhop)
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
outf = 5 # number of iterations between output files
outf = 5
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Outer",0,lambda x : -x[:, 1]),("Outer",1,lambda x : x[:, 0]),("Inner",0,0.),("Inner",1,0.),("PtFix",2,0.)]
#g is define in fluid_problem
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)
tic = time.clock()
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
#Computation of the new velocities
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
#number of sub time step
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
#Contact solver
for i in range(nsub) :
p.iterate(dt/nsub, forces,1e-6)
tol = 1e-6
p.iterate(dt/nsub, forces, tol)
#Localisation of the grains in the fluid
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
t += dt
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
......
testcases/depot-2d/dep.py
View file @ fdce8d90
......@@ -8,9 +8,15 @@ import time
import shutil
import random
#grains creation and placement + physical boundaries definition
'''r is the radius of the grains
ly is the height of the box
rhop is the particles density'''
def genInitialPosition(filename, r, ly, rhop) :
p = scontact2.ParticleProblem()
#Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom","BottomOut","TopOut"])
#Definition of the points where the grains are located
x = np.arange(-0.025+r, 0.025-r, r)
y = np.arange(0.06-r, .04, -r)
x, y = np.meshgrid(x, y)
......@@ -18,6 +24,7 @@ def genInitialPosition(filename, r, ly, rhop) :
y = y.flat
for i in range(len(x)) :
rhop = random.choice([1100,1200,1300])
#Addition of an particle object at each point
p.add_particle((x[i], y[i]), r/2, r**2/4 * np.pi * rhop);
p.write_vtk(filename,0,0)
......@@ -29,62 +36,65 @@ if not os.path.isdir(outputdir) :
t = 0
ii = 0
r = 2e-4
ly = 5e-2
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters
g = -9.81
rho = 1000
rhop = 1500
nu = 1e-6
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 100
r = 2e-4 # radius
ly = 5e-2 # box height
g = -9.81 # gravity
rho = 1000 # fluid density
rhop = 1500 # grains density
nu = 1e-6 # kinematic viscosity
tEnd = 100 # final time
#numerical parameters
lcmin = 0.001 # approx r*100 but should match the mesh size
dt = 1e-2
alpha = 2.5e-5
epsilon = alpha*lcmin**2 /nu
lcmin = 0.001 # mesh size
dt = 1e-2 # time step
alpha = 2.5e-5 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parameter
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#Object particles creation
genInitialPosition(outputdir, r, ly, rhop)
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 20
outf1 = 100000
outf = 20 # number of iterations between output files
ii = 0
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Top",2,0.),("TopOut",1,0),("Top",1,0),("BottomOut",1,0),("Bottom",1,0),("Lateral",0,0.),("Lateral",1,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)
tic = time.clock()
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
#Computation of the new velocities
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
#number of sub time step
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
#Contact solver
for i in range(nsub) :
p.iterate(dt/nsub, forces,1e-6)
tol = 1e-6
p.iterate(dt/nsub, forces, tol)
#Localisation of the grains in the fluid
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
t += dt
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
......
testcases/depot-3d/depot.py
View file @ fdce8d90
......@@ -6,24 +6,26 @@ import numpy as np
import os
import time
import shutil
import random
#grains creation and placement + physical boundaries definition
#The radii of the grains are defined in radius.csv and are computed to fit the distribution size given in the article of Andre et al. (2011) "Real-time analysis of the growth of granular media by an ultrasonic method: Application to the sedimentation of glass balls in water"
def genInitialPosition(filename, rhop) :
p = scontact3.ParticleProblem()
#Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
#Reading the file containing the radius
myRadius = np.genfromtxt('radius.csv', delimiter=',')
r1=np.amin(myRadius)*1e-6
r2=np.amax(myRadius)*1e-6
i = 0
#Positioning the grains on a regular grid
for y in np.arange(0, -.004, -2*r2):
for z in np.arange(.01675-r2, -.01675, -2*r2):
for x in np.arange(0.014-r2, -.014+r2, -2*r2):
p.add_particle((x, y, z), myRadius[i%500000]*1e-6, 4./3.* (myRadius[i%500000]*1e-6)**3 * np.pi * rhop);
i += 1
#Translation of the grains to the top of the box
p.position()[:, 1] += .15/2
print(len(p.position()[:, 1]))
print("Number of grains=%d" % len(p.position()[:, 1]))
p.write_vtk(filename, 0, 0)
......@@ -35,58 +37,57 @@ t = 0
ii = 0
p = scontact3.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters
g = -9.81
rho = 1000
rhop = 2640
nu = 1e-6
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 100
g = -9.81 # gravity
rho = 1000 # fluid density
rhop = 2640 # grains density
nu = 1e-6 # kinematic viscosity
tEnd = 100 # final time
#numerical parameters
lcmin = 0.001 # approx r*100 but should match the mesh size
dt = 1e-3
alpha = 2.5e-4
epsilon = alpha*lcmin**2 /nu#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5
lcmin = 0.001 # mesh size
dt = 1e-3 # time step
alpha = 2.5e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
#Object particles creation
genInitialPosition(outputdir, rhop)
p = scontact3.ParticleProblem()
p.read_vtk(outputdir,0)
#print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 100
outf = 100 # number of iterations between output files
strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]#("Top",3,0.),("Top",1,0.),("Top",0,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.)]
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal-width velocity; field 1 is vertical velocity; field 2 is the horizontal-depth velocity; field 3 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
print("Initialisation fluid done!\n")
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
print("Initialisation particles done!\n")
tic = time.clock()
fluid.export_vtk(outputdir,0,0)
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
#Computation of the new velocities
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
#number of sub time step
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
#Contact solver
for i in range(nsub) :
p.iterate(dt/nsub, forces)
tol = 1e-6
p.iterate(dt/nsub, forces, tol)
#Localisation of the grains in the fluid
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
t += dt
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
......
testcases/drop-2d/InteractingDrops/Diag/2DiagDrops.py
View file @ fdce8d90
......@@ -8,24 +8,35 @@ import time
import shutil
import random
#Physical parameters for the drops are the ones presented by Machu et al. (2001) "Coalescence, torus formation and breakup of sedimenting drops: experiments and computer simulations"
#grains creation and placement + physical boundaries definition
'''filename is the name of the output file
r is the radius of the grains
rout is the radius of the drop
rhop is the particles density
compacity is the solid volume fraction inside the drop'''
def genInitialPosition(filename, r, rout, rhop, compacity) :
p = scontact2.ParticleProblem()
#Loading of the mesh.msh file specifying physical boundaries name
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","Lateral"])
#Space between the grains to obtain the expected compacity
N = compacity*4*rout**2/(np.pi*r**2)
e = 2*rout/(N)**(1./2.)
#Definition of the points where the grains are located
x = np.arange(rout, -rout, -e)
x, y = np.meshgrid(x, x)
R2 = x**2 + y**2
keep = R2 < (rout - r)**2
x = x[keep]
y = y[keep]
#First drop
for i in range(x.shape[0]) :
p.add_particle((x[i]+random.uniform(-e/2+r,e/2-r)+rout, y[i]+random.uniform(-e/2+r,e/2-r)-2*rout), r, r**2 * np.pi * rhop);
#Second drop
for i in range(x.shape[0]) :
p.add_particle((x[i]+random.uniform(-e/2+r,e/2-r)-rout, y[i]+random.uniform(-e/2+r,e/2-r)+2*rout), r, r**2 * np.pi * rhop);
#Vertical shift of the grains to the top of the box
p.position()[:, 1] += 0.22
p.write_vtk(filename,0,0)
......@@ -37,69 +48,66 @@ if not os.path.isdir(outputdir) :
t = 0
ii = 0
#Enable computation of the initial velocity
initialize = 1
r=25e-6
R = 2.7e-3
compacity = 0.03
drho = 35.4
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
#physical parameters
g = -9.81
rhop = 2400
nu = 1e-4
rho = rhop-drho/compacity
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 1000
g = -9.81 # gravity
rhop = 2400 # grains density
nu = 1e-4 # kinematic viscosity
drho = 35.4 # density differences rhop-rho
compacity = 0.03 # solid volume fraction in the drop
rho = rhop-drho/compacity # fluid density
r=25e-6 # grains radii
R = 2.7e-3 # drop radius
mu = nu*rho # dynamic viscosity
tEnd = 1000 # final time
print("RHOP = %g" % rhop)
#numerical parameters
lcmin = 0.0008 # approx r*100 but should match the mesh size
dt = 5e-2
alpha = 1e-3
epsilon = alpha*lcmin**2 /nu
lcmin = 0.0008 # mesh size
dt = 5e-2 # time step
alpha = 1e-3 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
mu = nu*rho
#Object particles creation
genInitialPosition(outputdir, r, R, rhop,compacity)
p = scontact2.ParticleProblem()
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 1
outf1 = 100000
outf = 1 # number of iterations between output files
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)