Skip to content

GitLab

Commit 2b51e93d authored by Nathan Coppin's avatar Nathan Coppin
Browse files

preparing release

parent d9667408
Pipeline #7810 passed with stages
in 3 minutes and 13 seconds
Hide whitespace changes
Inline Side-by-side
testcases/couette-2d/couette.py
View file @ 2b51e93d
......@@ -104,7 +104,6 @@ else :
p.read_vtk(outputdir,0)
p.set_friction_coefficient(0.1,"Sand","Sand") # Particle-Particle
p.set_friction_coefficient(0.1,"Sand","Steel")# Particle-Wall
p.set_tangent_boundary_velocity("Inner",-0.1)
# Initial time and iteration
t = 0
ii = 0
......
testcases/drag-2d/depot/depot.py
View file @ 2b51e93d
......@@ -55,12 +55,12 @@ def genInitialPosition(filename, r, H, ly, lx, rin, rhop) :
p.write_vtk(filename,0,0)
# Define output directory
outputdir = "outputPoly"
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = np.array([0,-9.81]) # gravity
r = 3e-3 # particles radius
r = 8e-3 # particles radius
rin = 0.01905 # central disk radius
rhop = 2500 # particles density
# Numerical parameters
......@@ -68,29 +68,26 @@ outf = 1000 # number of iterations between out
dt = 1e-3 # time step
tEnd = 10 # final time
# Geometrical parameters
ly = 0.7 # particles area height
lx = 0.15 # particles area widht
H = 2 # domain height
ly = 0.7 # particles area height
lx = 0.15 # particles area widht
H = 2 # domain height
#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, H, ly, lx, rin, rhop)
p = scontact.ParticleProblem(2,True,True)
p.set_friction_coefficient(0.15,"Sand","Sand")
p.set_friction_coefficient(0.1,"Sand","Steel")
p.set_friction_coefficient(0.3,"Sand","Sand")
p.set_friction_coefficient(0.3,"Sand","Steel")
p.read_vtk(outputdir,0)
print("Number of particles : %d"%p.n_particles())
# Initial time and iteration
t = 0
ii = 0
G = np.zeros_like(p.velocity())
G[:,1] = p.mass()[:,0]*g[1]
#
# COMPUTATION LOOP
#
while t < tEnd :
time_integration.iterate(None, p, dt, min_nsub=1,contact_tol = 5e-7,external_particles_forces=G)
time_integration.iterate(None, p, dt, min_nsub=1,contact_tol = 5e-7,external_particles_forces=g*p.mass())
t += dt
# Output files writting
if ii %outf == 0 :
......
testcases/drag-2d/drag.py
View file @ 2b51e93d
......@@ -32,12 +32,11 @@ import time
import shutil
import random
# Physical parameters
mupart = 0.
vit = -0.336 #stream velocity
vit = -0.1 #stream velocity
rhop = 2500 #particles density
rho = 1000 #fluid density
nu = 1e-6 #fluid kinematic viscosity
g = np.array([0,-9.81]) #gravity
g = np.array([0,-9.81]) #gravity
# Numerical parameters
dt = 1e-3 #time step
t = 0 #initial time
......@@ -45,38 +44,25 @@ tEnd = 13 #final time
ii = 0 #iteration number
outf = 10 #iterations between data frames
# Define output directory
outputdir = "output" + str(abs(vit)).replace(".","_") + str(mupart).replace(".","_")
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
filename = outputdir + "/Drag.csv"
if os.path.exists(filename):
print("Remove old results? (y/n)\n")
decision = input()
if decision.lower() == 'y':
os.remove(filename)
#
# PARTICLE PROBLEM
#
#Injected particles
p = scontact.ParticleProblem(2,True,True)
p.read_vtk("depot/output2",11)
p.read_vtk("depot/output",10)
p.remove_particles_flag((p.position()[:,1]-p.r()[:,0]>0.5)*(p.position()[:,1]+p.r()[:,0] <0.53))
p.contact_forces()[:] = 0
p.velocity()[:] = 0
p.omega()[:] = 0
#Flow particles
p2 = scontact.ParticleProblem(2,True,True)
p2.read_vtk("depot/output2",11)
p2.read_vtk("depot/output",10)
p2.remove_particles_flag(p2.position()[:,1]+p2.r()[:,0] <0.5)
p2.contact_forces()[:] = 0
p2.velocity()[:] = 0
p2.omega()[:] = 0
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner"],material = "Steel")
p2.load_msh_boundaries("mesh.msh", ["Left", "Right"],material = "Plexi")
p2.set_friction_coefficient(mupart,"Sand","Sand")
p2.set_friction_coefficient(0.,"Sand","Steel")
p2.set_friction_coefficient(0.,"Sand","Plexi")
p2.load_msh_boundaries("mesh.msh", ["Inner","Left","Right"],material = "Steel")
p2.set_friction_coefficient(0.3,"Sand","Sand")
p2.set_friction_coefficient(0.3,"Sand","Steel")
p2.write_vtk(outputdir, 0, 0)
#
# FLUID PROBLEM
......@@ -94,29 +80,25 @@ fluid.export_vtk(outputdir, 0.,0)
# COMPUTATION LOOP
#
def accumulate(F):
F+=np.sum(p2.get_boundary_impulsion("Inner"),axis=0)
F+=np.sum(p2.get_boundary_forces("Inner"),axis=0)
while t < tEnd :
p2.forced_flag()[(p2.position()[:,1] + p2.r()[:,0] <= -0.5)] = 1
if t > 3 or t<10*dt:
p2.velocity()[p2.forced_flag()==1,:] = [0, vit]
else :
p2.velocity()[p2.forced_flag()==1,:] = [0, 0]
p2.velocity()[p2.forced_flag()==1,:] = [0, vit]
p2.omega()[p2.forced_flag()==1] = 0
if max(p2.position()[:,1]+p2.r()[:,0]) < 0.5 and t < 25 :
p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity()[:,:] + [0,vit],tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity()[:,:] + [0,vit],tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.55))
F = np.zeros(2)
#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 subdt : accumulate(F))
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))
with open(filename,"a") as file1:
F /= dt
#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-2d/mesh.geo
View file @ 2b51e93d
rin = 0.01905;
wout = 0.56;
lout = 0.15;
lcin = wout/20;
lcout = wout/20;
lcin = wout/2;
lcout = wout/2;
Point(1) = {0, 0, 0};
......
testcases/drag-3d/depot/depot.py
View file @ 2b51e93d
......@@ -35,7 +35,6 @@ import time
import shutil
import random
import sys
numDep = int(sys.argv[1])
def genInitialPosition(filename, r, H, ly, lx, lz, rin, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure
......@@ -64,30 +63,30 @@ def genInitialPosition(filename, r, H, ly, lx, lz, rin, rhop) :
p.add_particle((x[i], y[j], z[k]), r1, (4/3)*r1**3 * np.pi * rhop,"Sand")
p.write_vtk(filename,0,0)
# Define output directory
outputdir = str(numDep)
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = np.array([0,-9.81,0]) #gravity
r = 3e-3 #particles radius
r = 8e-3 #particles radius
rin = 0.01905 #central cylinder radius
rhop = 2500 #particles density
# Numerical parameters
outf = 5000 #number of iterations between output files
outf = 50 #number of iterations between output files
dt = 1e-3 #time step
tEnd = 5 + 5*dt #final time
tEnd = 10 #final time
# Geometrical parameters
ly = 0.7 #particles area height
lx = 0.15 #particles area widht
ly = 0.7 #particles area height
lx = 0.15 #particles area widht
lz = 0.05 #particles area depth
H = 1.4 #domain height
H = 1.3 #domain height
#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, H, ly, lx, lz, rin, rhop)
p = scontact.ParticleProblem(3,True,True)
p.set_friction_coefficient(0.15,"Sand","Sand")
p.set_friction_coefficient(0.3,"Sand","Sand")
p.set_friction_coefficient(0.3,"Sand","Steel")
p.read_vtk(outputdir,0)
# Initial time and iteration
......
testcases/drag-3d/drag.py
View file @ 2b51e93d
......@@ -32,44 +32,35 @@ import time
import shutil
import random
import sys
listeVit = [14.39,36.32,64.06,95.57,130.75,169.32,211.05,255.63,303.22,354.05,408.80,467.48,530.59]
inter = int(sys.argv[1])
index = int(inter/5)
numDep = inter%5
# Physical parameters
vit = -0.001*listeVit[index] #stream velocity
vit = -0.05 #stream velocity
muwall = 0.3 #friction coefficient between the walls and the particles
mupart = 0.15 #friction coefficient between particles
rho = 1000 #fluid density
nu = 5e-5 #fluid kinematic viscosity
nu = 1e-6 #fluid kinematic viscosity
g = np.array([0,-9.81,0]) #gravity
# Numerical parameters
dt = 1e-3 #time step
t = 0 #initial time
tEnd = 10 #final time
ii = 0 #iteration number
outf = 100 #iterations between data frames
outf = 100 #iterations between data frames
# Define output directory
outputdir = "drag" + str(int(listeVit[index])) + "/" + str(numDep)
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
filename = outputdir + "/Drag.csv"
if os.path.exists(filename):
print("Remove old results? (y/n)\n")
decision = input()
if decision.lower() == 'y':
os.remove(filename)
#
# PARTICLE PROBLEM
#
#Injected particles
p = scontact.ParticleProblem(3,True,True)
p.read_vtk("outputhahalol",5)
p.read_vtk("depot/output",10)
p.remove_particles_flag((p.position()[:,1]-p.r()[:,0]>0.5)*(p.position()[:,1]+p.r()[:,0] <0.52))
p.velocity()[:,1] = vit
#Flow particles
p2 = scontact.ParticleProblem(3,True,True)
p2.read_vtk("outputhahalol",5)
p2.read_vtk("depot/output",10)
p2.remove_particles_flag(p2.position()[:,1]+p2.r()[:,0] <0.5)
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner","Left", "Right","Front","Rear"],material = "Steel")
......@@ -94,7 +85,7 @@ fluid.export_vtk(outputdir, 0.,0)
# COMPUTATION LOOP
#
def accumulate(F):
F+=np.sum(p2.get_boundary_impulsion("Inner"),axis=0)
F+=np.sum(p2.get_boundary_forces("Inner"),axis=0)
while t < tEnd :
p2.forced_flag()[(p2.position()[:,1] + p2.r()[:,0] <= -0.5)] = 1
p2.velocity()[p2.forced_flag()==1,:] = [0, vit, 0]
......@@ -107,8 +98,7 @@ while t < tEnd :
time_integration.iterate(fluid,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda subdt : accumulate(F))
with open(filename,"a") as file1:
F += fluid.boundary_force()[0,:]
meanVel = abs((np.mean(np.linalg.norm(p2.velocity()[(p2.position()[:,1] > -0.45)*(p2.position()[:,1]<-0.3),:2],axis=1)) + np.mean(np.linalg.norm(p2.velocity()[(p2.position()[:,1] < 0.45)*(p2.position()[:,1]>0.3),:2],axis=1)))/2)
file1.write( str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+";"+str(meanVel)+"\n")
file1.write( str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+"\n")
t += dt
# Output files writing
if ii %outf == 0 :
......
testcases/hopper-2d/hopper.py
View file @ 2b51e93d
......@@ -26,7 +26,6 @@ from migflow import fluid
from migflow import time_integration
import numpy as np
import os
num = "05"
def genInitialPosition(filename, r, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure.
......@@ -40,9 +39,9 @@ def genInitialPosition(filename, r, rhop) :
#Particles structure builder
p = scontact.ParticleProblem(2, friction_enabled = True,rotation_enabled=True)
#Load mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh" + num + ".msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh" + num + ".msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"], material = "Plexi")
p.load_msh_boundaries("mesh" + num + ".msh", ["LockDown"], material = "Steel")
p.load_msh_boundaries("mesh.msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh.msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"], material = "Plexi")
p.load_msh_boundaries("mesh.msh", ["LockDown"], material = "Steel")
#Particles centre are placed on a regular grid
# Add a grain at each centre position
for y in np.arange(1+r, 2, 2*r*1.05):
......@@ -52,7 +51,7 @@ def genInitialPosition(filename, r, rhop) :
p.write_vtk(outputdir, 0, 0)
# Physical parameters
g = np.array([0,-9.81])
r = 3e-3
r = 8e-3
rhop = 2500
mupart= 0.12
# Numerical parameters
......@@ -60,18 +59,13 @@ dt = 1e-3 #time step
t = 0 #initial time
tEnd = 13 #final time
ii = 0 #iteration number
outf = 1000 #iterations between data frames
outf = 1000 #iterations between data frames
# Define output directory
outputdir = "output" + num + str(mupart).replace(".","_")
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#Creating file for drag data
filename = outputdir + "/Drag.csv"
if os.path.exists(filename):
print("Remove old results? (y/n)\n")
decision = input()
if decision.lower() == 'y':
os.remove(filename)
#
# PARTICLE PROBLEM
#
......@@ -86,9 +80,9 @@ p2.remove_particles_flag((p2.position()[:,0] + p2.r()[:,0] <0.15))
p2.remove_particles_flag((p2.position()[:,0] + p2.r()[:,0] >-0.15))
p2.position()[:,1] -= 0.6
p.clear_boundaries()
p.load_msh_boundaries("mesh" + num + ".msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh" + num + ".msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"], material = "Plexi")
p.load_msh_boundaries("mesh" + num + ".msh", ["LockDown"], material = "Steel")
p.load_msh_boundaries("mesh.msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh.msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"], material = "Plexi")
p.load_msh_boundaries("mesh.msh", ["LockDown"], material = "Steel")
p.set_friction_coefficient(mupart,"Sand","Sand")
p.set_friction_coefficient(0.1,"Sand","Steel")
p.set_friction_coefficient(0.1,"Sand","Plexi")
......@@ -98,13 +92,13 @@ p.write_vtk(outputdir, 0, 0)
#
opened = 0
def accumulate(F):
F+=np.sum(p.get_boundary_impulsion("Inner"),axis=0)
F+=np.sum(p.get_boundary_forces("Inner"),axis=0)
while t < tEnd :
if t>=3 and opened == 0:
opened = 1
p.clear_boundaries()
p.load_msh_boundaries("mesh" + num + ".msh", ["Inner"],material="Steel")
p.load_msh_boundaries("mesh" + num + ".msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"],material="Plexi")
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(p.position()[:,1])+max(p.r()) <0.5) and t > 3 and t < 25:
p.add_particles(p2.position(),p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
......@@ -114,12 +108,9 @@ while t < tEnd :
#Removing particles outside the hopper
p.remove_particles_flag( (p.position()[:,1] >-0.6))
#Computing mean velocity and writing to file
speedAbove = np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1] > 0.3),1])
speedBelow = np.mean(p.velocity()[(p.position()[:,1] > -0.45)*(p.position()[:,1] < -0.3),1])
with open(filename,"a") as file1:
file1.write(str(F[0])+";"+str(F[1])+";"
+str(t)+";"+str(0 if np.isnan(speedAbove) else speedAbove)+";"
+str(0 if np.isnan(speedBelow) else speedBelow)+"\n")
+str(t)+"\n")
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
......
testcases/hopper-2d/hopperFluid.py deleted 100644 → 0
View file @ d9667408
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
# TESTCASE DESCRIPTION
# Study of the drag on a fixed sphere in a immersed granular flow
from migflow import scontact
from migflow import fluid
from migflow import time_integration
import numpy as np
import os
import time
import shutil
import random
numdep = "004"
# Physical parameters
rapport = 2
g = np.array([0,-9.81])
rho = 1000 #fluid density
nu = 1e-6 #fluid kinematic viscosity
# Numerical parameters
dt = 1e-3 #time step
t = 0 #initial time
tEnd = 10 #final time
ii = 0 #iteration number
outf = 100 #iterations between data frames
# Define output directory
inputdir = "depot/depot" + numdep
outputdir = "water-" + numdep + "-" + str(rapport).replace(".","_")
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#Creating file for drag data
filename = outputdir + "/Drag.csv"
if os.path.exists(filename):
print("Remove old results? (y/n)\n")
decision = input()
if decision.lower() == 'y':
os.remove(filename)
#
# PARTICLE PROBLEM
#
p = scontact.ParticleProblem(2,True, True)
p.read_vtk(inputdir,10)
p2 = scontact.ParticleProblem(2,True, True)
p2.read_vtk(inputdir,10)
p2.contact_forces()[:] = 0
p2.omega()[:] = 0
p2.velocity()[:,:] = 0
p2.mass()[:] *=rapport/2.5
p.contact_forces()[:] = 0
p.omega()[:] = 0
p.velocity()[:,:] = 0
p.mass()[:] *=rapport/2.5
p2.remove_particles_flag((p2.position()[:,1] + p2.r()[:,0] <1.2))
p2.remove_particles_flag((p2.position()[:,1] - p2.r()[:,0] >1.15))
p2.position()[:,1] += 0.05
flipped = np.copy(p2.position())
flipped[:,0] = -flipped[:,0]
p.clear_boundaries()
p.load_msh_boundaries(inputdir+"/mesh.msh", ["Inner"],material="Steel")
p.load_msh_boundaries(inputdir+"/mesh.msh", ["RightUp","Right","RightDown","LeftDown","Left","LeftUp"],material="Plexi")
p.load_msh_boundaries(inputdir+"/mesh.msh", ["LockUp","LockDown"],material="Steel")
p.remove_particles_flag((p.position()[:,1] + p.r()[:,0] <1.2))
p.add_particles(p2.position(),p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
p.add_particles(p2.position()[:,:]+[0,0.05],p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
p.add_particles(p2.position()[:,:]+[0,0.1],p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
p.add_particles(p2.position()[:,:]+[0,0.15],p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
p.add_particles(p2.position()[:,:]+[0,0.2],p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
p.set_friction_coefficient(0.3,"Sand","Sand")
p.set_friction_coefficient(0.12,"Sand","Steel")
p.set_friction_coefficient(0.05,"Sand","Plexi")
p.write_vtk(outputdir, 0, 0)
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=0)
#Set the mesh geometry for the fluid computation
fluid.load_msh(inputdir+"/mesh.msh")
fluid.set_wall_boundary("Inner")
fluid.set_wall_boundary("Left")
fluid.set_wall_boundary("LeftDown")
fluid.set_wall_boundary("LeftUp")
fluid.set_wall_boundary("Right")
fluid.set_wall_boundary("RightDown")
fluid.set_wall_boundary("RightUp")
fluid.set_open_boundary("LockDown",pressure = -rho*g*2.1)
fluid.set_open_boundary("LockUp",pressure = 0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces(),init=True)
fluid.export_vtk(outputdir, 0.,0)
#
# COMPUTATION LOOP
#
def accumulate(F):
F+=np.sum(p.get_boundary_impulsion("Inner"),axis=0)
while t < tEnd :
#Adding new particles
if (max(p.position()[:,1])+max(p.r()) <1.4):
p.add_particles(p2.position()[:,:]+[0,0.2] if random.choice([True,False]) else flipped[:,:]+[0,0.2],p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
#Solving the contacts
F = np.zeros(2)
time_integration.particle_changed(fluid,p)
time_integration.iterate(fluid,p,dt,min_nsub=5,contact_tol=2.5e-7,external_particles_forces=g*p.mass(),after_sub_iter = lambda subdt : accumulate(F))
#Removing particles outside the hopper
p.remove_particles_flag( (p.position()[:,1] >-0.59))
#Computing mean velocity and writing to file
speedAbove = np.mean(p.velocity()[(p.position()[:,1] < 0.45)*(p.position()[:,1] > 0.3),1])
speedBelow = np.mean(p.velocity()[(p.position()[:,1] > -0.45)*(p.position()[:,1] < -0.3),1])
#print("Speed Above",speedAbove,"Speed Below",speedBelow)
with open(filename,"a") as file1:
F += fluid.boundary_force()[0,:];
file1.write(str(F[0])+";"+str(F[1])+";"+str(t)+";"
+str(0 if np.isnan(speedAbove) else speedAbove)+";"+str(0 if np.isnan(speedBelow) else speedBelow)+"\n")
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
testcases/hopper-3d/hopper.py
View file @ 2b51e93d
......@@ -26,24 +26,9 @@ from migflow import fluid
from migflow import time_integration
from migflow import lmgc90Interface
from pylmgc90 import pre
import numpy as np
import os
import sys
#liste = ["015","0175","02","0225","025","03","04","05","06","07","08","09","10","11","12","13","14"]
#inter = int(sys.argv[1])
#index = int(inter/5)
#numM = liste[index]
#num = inter%5
use_lmgc = False
num = int(sys.argv[1])
friction= [[0.0,0.0],[0.0,0.1],[0.0,0.2],[0.0,0.3],[0.0,0.4],[0.0,0.5],
[0.1,0.0],[0.1,0.1],[0.1,0.2],[0.1,0.3],[0.1,0.4],[0.1,0.5],
[0.2,0.0],[0.2,0.1],[0.2,0.2],[0.2,0.3],[0.2,0.4],[0.2,0.5],
[0.3,0.0],[0.3,0.1],[0.3,0.2],[0.3,0.3],[0.3,0.4],[0.3,0.5],
[0.4,0.0],[0.4,0.1],[0.4,0.2],[0.4,0.3],[0.4,0.4],[0.4,0.5],
[0.5,0.0],[0.5,0.1],[0.5,0.2],[0.5,0.3],[0.5,0.4],[0.5,0.5]]
def genInitialPosition(filename, r, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure.
......@@ -75,19 +60,13 @@ dt = 1e-3 #time step
t = 0 #initial time
tEnd = 40 #final time
ii = 0 #iteration number
outf = 100 #iterations between data frames
outf = 1000 #iterations between data frames
# Define output directory
#outputdir ="/CECI/home/ucl/mema/ncoppin/hopper" + numM + "/" + str(num)
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#Creating file for drag data
filename = outputdir + "/Drag.csv"
if os.path.exists(filename):
print("Remove old results? (y/n)\n")
decision = input()
if decision.lower() == 'y':
os.remove(filename)
#
# PARTICLE PROBLEM
#
......@@ -99,8 +78,8 @@ if use_lmgc :
p = lmgc90Interface.ParticleProblem(3)
else :
p = scontact.ParticleProblem(3,friction_enabled = True,rotation_enabled=True)
p.set_friction_coefficient(friction[num][0],"Sand","Sand")
p.set_friction_coefficient(friction[num][1],"Sand","Steel")
p.set_friction_coefficient(0.3,"Sand","Sand")
p.set_friction_coefficient(0.3,"Sand","Steel")
p2 = scontact.ParticleProblem(3,friction_enabled = True,rotation_enabled=True)
p.read_vtk(outputdir, 0)