Skip to content

GitLab

Commit 2be74e44 authored by Hugo Ginestet's avatar Hugo Ginestet
Browse files

add testcases

parent f0d9da13
Pipeline #5429 passed with stage
in 27 seconds
Hide whitespace changes
Inline Side-by-side
testcases/hugo/adaptation1.py 0 → 100644
View file @ 2be74e44
from migflow import fluid
from migflow import scontact
import copy
import numpy as np
import math
import os
import time
import shutil
import random
L = 10.
lc = 0.2
h_min =0.25
h_max= 0.25
r = 1.0
para = [L,lc,h_min,h_max,[0.,0.,r]]
###################################################################################################################
"""
Création de la géométrie pour la création du maillage
"""
def creation_point(fic,Point):
fic.write('Point('+str(Point[0])+') = {'+str(Point[1])+', '+str(Point[2])+', 0.0, '+str(Point[3])+'};\n')
def creation_ligne(fic,i,Point1,Point2):
fic.write('Line('+str(i)+') = {'+str(Point1)+', '+str(Point2)+'};\n')
def creation_loop(fic,i,Ligne1,Ligne2,Ligne3,Ligne4):
fic.write('Line Loop('+str(i)+') = {'+str(Ligne1)+', '+str(Ligne2)+', '+str(Ligne3)+', ' +str(Ligne4)+'};\n\n')
def initialisation(fic,h_min,h_max):
fic.write('SetFactory("OpenCASCADE");\n')
fic.write('Mesh.CharacteristicLengthMin ='+ str(h_min) +' ;\n')
fic.write('Mesh.CharacteristicLengthMax ='+ str(h_max) +' ;\n\n')
fic.write('Mesh.MinimumCirclePoints = 1 ;\n\n')
fic.write('Mesh.MinimumCurvePoints = 2 ;\n\n')
def creation_box(fic,L,lc):
Point = [[1,-L/2.,-L/2.,lc],[2,L/2.,-L/2.,lc],[3,L/2.,L/2.,lc],[4,-L/2.,L/2.,lc]]
creation_point(fic,Point[0])
for k in range(3):
creation_point(fic,Point[k+1])
creation_ligne(fic,k+1,Point[k][0],Point[k+1][0])
creation_ligne(fic,4,Point[-1][0],Point[0][0])
creation_loop(fic,1,Point[0][0],Point[1][0],Point[2][0],Point[3][0])
fic.write('Surface(1) = {1}; \n\n')
def creation_disk(fic,i,Point):
fic.write('Disk('+str(i)+') = {'+str(Point[0])+', '+str(Point[1])+', 0., '+str(Point[2])+'};\n\n')
def creation_surface(fic,nb):
fic.write('BooleanDifference{Surface{1};Delete;}{Surface{2};Delete;}\n\n')
def creation_contour_rectangle(fic):
fic.write('Physical Surface("Domain") = {1};\n')
fic.write('Physical Line("Bottom") = {6};\n')
fic.write('Physical Line("Top") = {9};\n')
fic.write('Physical Line("Lateral") = {7,8};\n')
fic.write('Physical Line("Grains") = {5};\n\n')
fic.write('Mesh.MshFileVersion = 2;')
def creation_maillage(name,para):
fic = open(name+'.geo',"w")
initialisation(fic,para[2],para[3])
creation_box(fic,para[0],para[1])
creation_disk(fic,2,para[4])
creation_surface(fic,1)
creation_contour_rectangle(fic)
fic.close()
os.system("gmsh -2 {}.geo".format(name))
#os.system("rm {}.geo".format(name))
creation_maillage('mesh',para)
para = [L,lc,h_min,h_max]
def depl_lineaire(pos,para,depl,pos_c,dt):
l = para[0]/2
#norm0 = math.sqrt((pos[0]-pos_c[0])**2 + (pos[1]-pos_c[1])**2)
velocity = [0, 0]
'''
if abs(norm0 - pos_c[2]) < 1e-4:
for i in range(2):
pos[i] += depl[i]
velocity[i] = depl[i]/dt
'''
for i in range(2):
if pos[i] <(pos_c[i] - pos_c[2]):
pos[i] += depl[i]*(pos[i]+l)/(pos_c[i]- pos_c[2]+l)
velocity[i] = (depl[i]*(pos[i]+l)/(pos_c[i]- pos_c[2]+l))/dt
elif pos[i] >(pos_c[i] + pos_c[2]):
pos[i] += depl[i]*(pos[i]-l)/(pos_c[i]+ pos_c[2]-l)
velocity[i] = (depl[i]*(pos[i]-l)/(pos_c[i]+ pos_c[2]-l))/dt
else:
pos[i]+= depl[i]
velocity[i] = depl[i]/dt
return velocity
def depl_lineaire2(pos,para,depl,pos_c):
l = para[0]/2
norm0 = math.sqrt((pos[0])**2 + (pos[1])**2)
if abs(pos[0]-pos_c[0]) > 1e-5:
x = (pos[1]-pos_c[1])/(pos[0]-pos_c[0])
pos0 = [norm0/math.sqrt(1+x**2),norm0*x/math.sqrt(1+x**2)]
else:
pos0 = [pos[0]-pos_c[0],pos[1]-pos_c[1]]
for i in range(2):
if pos[i] <pos_c[i]:
pos[i] += depl[i]*(pos0[i]+l)/(-pos_c[2]+l)
else:
pos[i] += depl[i]*(pos0[i]-l)/(pos_c[2]-l)
return
def calcul_norm(para,depl,pos_c):
l = para[0]/2
for k in range(2):
if pos_c[0] > 0:
if pos_c[1]>0:
return min(abs(pos_c[0]-l),abs(pos_c[1]-l))
else:
return min(abs(pos_c[0]-l),abs(pos_c[1]+l))
else:
if pos_c[1]>0:
return min(abs(pos_c[0]+l),abs(pos_c[1]-l))
else:
return min(abs(pos_c[0]+l),abs(pos_c[1]+l))
def depl_radial(pos,depl,pos_c,norm_max,para,dt):
r = pos_c[2]
norm0 = math.sqrt((pos[0]-pos_c[0])**2+(pos[1]-pos_c[1])**2)
if norm0 < norm_max:
for i in range(2):
pos[i] += depl[i]*(norm0 - norm_max)/(r - norm_max)
return [depl[0]*(norm0 - norm_max)/((r - norm_max)*dt), depl[1]*(norm0 - norm_max)/((r - norm_max)*dt)]
def depl_radial0(pos0,pos,depl,pos_c,norm_max,dt):
r = pos_c[2]
norm0 = math.sqrt((pos0[0]-pos_c[0])**2+(pos0[1]-pos_c[1])**2)
new_pos = copy.deepcopy(pos0)
velocity = [0,0]
if norm0 < norm_max:
for i in range(2):
new_pos[i] += depl[i]*(norm0 - norm_max)/(r - norm_max)
if abs(r-norm0)>1e-10:
velocity[i] = (new_pos[i] - pos[i])/dt
else:
for i in range(2):
new_pos[i] = pos[i]
return [velocity,new_pos]
def geo_ligne(nb,l,vec):
L = [[0,0]]
pas = l/nb
norm =math.sqrt(vec[0]**2+vec[1]**2)
for k in range(nb):
L.append([pas*vec[0]/norm,pas*vec[1]/norm])
return L
def geo_cercle(nb,r):
theta = 2*math.pi/nb
L = [[0,0]]
for k in range(1,nb):
toto = [r*math.cos(k*theta) - r*math.cos((k-1)*theta), r*math.sin(k*theta) -r*math.sin((k-1)*theta)]
L.append(toto)
return L
def retourne_zero(pos_c,nb):
return [[(0.-pos_c[0])/nb,(0.-pos_c[1])/nb] for i in range(nb)]
#########################################################################################
outputdir = "output0"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = -9.81 # gravity
# particles radius
rhop = 1500 # particles density
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
# Numerical parameters
outf = 2 # number of iterations between output files
dt = 0.5e-2 # time step
tEnd = 5
# final time
nbtot = int(tEnd/dt)
#PARTICLE PROBLEM
#
def genInitialPosition(filename,r,L,H,para,rhop):
# Particles structure builder
p = scontact.ParticleProblem(2)
# Load mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom","Grains"])
p.add_particle((5., 5.), r, r**2 * np.pi * rhop)
p.write_vtk(filename,0,0)
#p = scontact.ParticleProblem(2)
#p.read_vtk(outputdir,0)
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
#fluid.set_wall_boundary("Grains", velocity = [10., 0.])
#fluid.set_open_boundary("Grains",velocity=[0.,0.],physical = True)
fluid.set_strong_boundary("Grains",0,10.)
fluid.set_strong_boundary("Grains",1,0.)
fluid.set_wall_boundary("Top",pressure=0)
fluid.export_vtk(outputdir,0,0)
tic = time.time()
rayon = 2.
nb1 = int(nbtot/(2*(np.pi+2)))
nb2 = nbtot -2*nb1
print((nb1,nb2))
liste_depl = geo_ligne(nb1,rayon,[1,0]) + 1*geo_cercle(nb2,rayon) + geo_ligne(nb1,rayon,[-1,0])
#liste_depl = 1*(geo_ligne(nb2,rayon,[0,-5])+geo_ligne(200,rayon,[2,1]))
#liste_depl = geo_ligne(nb1,rayon,[0,-1])
n = len(liste_depl)
N = len(fluid.coordinates())
pos_centre = [0.,0.,r]
depl = copy.deepcopy(pos_centre)
#
# COMPUTATION LOOP
#
k = 0
t = 0
ii = 0
'''
pos0= copy.deepcopy(fluid.coordinates())
cl_vitesse = [0.,0.]
for k in range(n):
#norm_max = calcul_norm(para,liste_depl[k],pos_centre)
for j in range(N):
#[poubelle,fluid.coordinates()[j]] = depl_radial0(pos0[j],fluid.coordinates()[j],depl,pos_centre,L/2,dt)
[fluid.velocity_mesh()[j],fluid.coordinates()[j]] = depl_radial0(pos0[j],fluid.coordinates()[j],depl,pos_centre,L/2,dt)
cl_vitesse[0] = liste_depl[k][0]/dt
cl_vitesse[1] = liste_depl[k][1]/dt
print(cl_vitesse)
#fluid.set_open_boundary("Grains", velocity=cl_vitesse)
fluid.set_strong_boundary("Grains",0,cl_vitesse[0])
fluid.set_strong_boundary("Grains",1,cl_vitesse[1])
fluid.implicit_euler(dt,newton_max_it=10)
fluid.export_vtk(outputdir, (k+1), (k+1))
depl[0] += liste_depl[k][0]
depl[1] += liste_depl[k][1]
print('{}\n'.format(k))
cl_vitesse[0] = 0
cl_vitesse[1] = 0
fluid.set_open_boundary("Grains", velocity=cl_vitesse)
'''
for k in range(n):
fluid.implicit_euler(dt,newton_max_it=10)
fluid.export_vtk(outputdir, (n+k+1), (n+k+1))
print('{}\n'.format(k+n))
'''
liste_depl = retourne_zero(pos_centre,30)
nb = len(liste_depl)
for k in range(nb):
norm_max = calcul_norm(para,liste_depl[k],pos_centre)
for j in range(N):
#fluid.velocity_mesh()[j] = depl_radial(fluid.coordinates()[j],liste_depl[k],pos_centre,norm_max,para,dt)
fluid.velocity_mesh()[j] = depl_lineaire(fluid.coordinates()[j],para,liste_depl[k],pos_centre,dt)
#depl_lineaire2(fluid.coordinates()[j],para,liste_depl[k],pos_centre)
#fluid.implicit_euler(dt)
fluid.export_vtk(outputdir, (n+k+1), (n+k+1))
pos_centre[0] += liste_depl[k][0]
pos_centre[1] += liste_depl[k][1]
print('{}\n'.format(k))
os.system("paraview {}/fluid.pvd".format(outputdir))
'''
\ No newline at end of file
testcases/hugo/chute.py 0 → 100644
View file @ 2be74e44
from migflow import fluid
from migflow import scontact
from migflow import fluid
from migflow import scontact
import copy
import numpy as np
import math
import os
import time
import shutil
import random
L = 5.
lc = 0.2
h_min = 0.25
h_max=0.25
r = 0.5
para = [L,lc,h_min,h_max,[0.,0.,r]]
###################################################################################################################
"""
Création de la géométrie pour la création du maillage
"""
def creation_point(fic,Point):
fic.write('Point('+str(Point[0])+') = {'+str(Point[1])+', '+str(Point[2])+', 0.0, '+str(Point[3])+'};\n')
def creation_ligne(fic,i,Point1,Point2):
fic.write('Line('+str(i)+') = {'+str(Point1)+', '+str(Point2)+'};\n')
def creation_loop(fic,i,Ligne1,Ligne2,Ligne3,Ligne4):
fic.write('Line Loop('+str(i)+') = {'+str(Ligne1)+', '+str(Ligne2)+', '+str(Ligne3)+', ' +str(Ligne4)+'};\n\n')
def initialisation(fic,h_min,h_max):
fic.write('SetFactory("OpenCASCADE");\n')
fic.write('Mesh.CharacteristicLengthMin ='+ str(h_min) +' ;\n')
fic.write('Mesh.CharacteristicLengthMax ='+ str(h_max) +' ;\n\n')
fic.write('Mesh.MinimumCirclePoints = 1 ;\n\n')
fic.write('Mesh.MinimumCurvePoints = 2 ;\n\n')
def creation_box(fic,L,lc):
Point = [[1,-L/2.,-L/2.,lc],[2,L/2.,-L/2.,lc],[3,L/2.,L/2.,lc],[4,-L/2.,L/2.,lc]]
creation_point(fic,Point[0])
for k in range(3):
creation_point(fic,Point[k+1])
creation_ligne(fic,k+1,Point[k][0],Point[k+1][0])
creation_ligne(fic,4,Point[-1][0],Point[0][0])
creation_loop(fic,1,Point[0][0],Point[1][0],Point[2][0],Point[3][0])
fic.write('Surface(1) = {1}; \n\n')
def creation_disk(fic,i,Point):
fic.write('Disk('+str(i)+') = {'+str(Point[0])+', '+str(Point[1])+', 0., '+str(Point[2])+'};\n\n')
def creation_surface(fic,nb):
fic.write('BooleanDifference{Surface{1};Delete;}{Surface{2};Delete;}\n\n')
def creation_contour_rectangle(fic):
fic.write('Physical Surface("Domain") = {1};\n')
fic.write('Physical Line("Bottom") = {6};\n')
fic.write('Physical Line("Top") = {9};\n')
fic.write('Physical Line("Lateral") = {7,8};\n')
fic.write('Physical Line("Grains") = {5};\n\n')
fic.write('Mesh.MshFileVersion = 2;')
def creation_maillage(name,para):
fic = open(name+'.geo',"w")
initialisation(fic,para[2],para[3])
creation_box(fic,para[0],para[1])
creation_disk(fic,2,para[4])
creation_surface(fic,1)
creation_contour_rectangle(fic)
fic.close()
os.system("gmsh -2 {}.geo".format(name))
#os.system("rm {}.geo".format(name))
creation_maillage('mesh',para)
para = [L,lc,h_min,h_max]
def depl_radial0(pos0,pos,depl,pos_c,norm_max,dt):
r = pos_c[2]
norm0 = math.sqrt((pos0[0]-pos_c[0])**2+(pos0[1]-pos_c[1])**2)
new_pos = copy.deepcopy(pos0)
velocity = [0,0]
if norm0 < norm_max:
for i in range(2):
new_pos[i] += depl[i]*(norm0 - norm_max)/(r - norm_max)
if abs(r-norm0)>1e-10:
velocity[i] = (new_pos[i] - pos[i])/dt
else:
for i in range(2):
new_pos[i] = pos[i]
return [velocity,new_pos]
def geo_ligne(nb,l,vec):
L = [[0,0]]
pas = l/nb
norm =math.sqrt(vec[0]**2+vec[1]**2)
for k in range(nb):
L.append([pas*vec[0]/norm,pas*vec[1]/norm])
return L
def genInitialPosition(filename, r, L, ly, lx, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure
Keyword arguments:
filename -- name of the output file
r -- max radius of the particles
H -- domain height
ly - particles area height
lx -- particles area width
rhop -- particles density
"""
# Particles structure builder
p = scontact.ParticleProblem(2)
# Load mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom","Grains"])
#Definition of the points where the particles are located
x = np.arange(-L/2+r, L/2-r, 2*r)
y = np.arange(-L/2+r, -L/2+ly-r, 2*r)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
# Add a grain at each centre position
for i in range(len(x)) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
p.write_vtk(filename,0,0)
# Define output directory
outputdir = "output2"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = -9.81 # gravity
rbille = 15.e-3 # particles radius
rhop = 1500 # particles
rho = 1000 # fluid density
nu = 1.e-6 # kinematic
# Numerical parameters
outf = 5 # number of iterations between output files
dt = 5.e-4 # time step
tEnd = 10 # final time
# Geometrical parameters
ly = L/3 # particles area height
lx = 4e-1 # particles area widht
# domain height
#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, rbille, L, ly, lx, rhop)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
# Initial time and iteration
t = 0
ii = 0
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
#fluid.set_wall_boundary("Grains")
fluid.set_wall_boundary("Top",pressure=0)
# Set location of the particles 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())
fluid.export_vtk(outputdir,0,0)
tic = time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
#
# COMPUTATION LOOP
#
rayon = 1.
nb1 = int(tEnd/dt*0.01)
print(nb1)
liste_depl = geo_ligne(nb1,rayon,[0.,-1])
n = len(liste_depl)
N = len(fluid.coordinates())
print(n)
pos_centre = [0.,0.,r]
depl = [pos_centre[0],pos_centre[1]]
pos0= copy.deepcopy(fluid.coordinates())
poubelle = [0., 0.]
k = 0
vitesse = [0.,0.]
v_part = [0., 0.]
while t < tEnd :
# Fluid solver
if k < (n-1):
#norm_max = calcul_norm(para,liste_depl[k],pos_centre)
for j in range(N):
[fluid.velocity_mesh()[j],fluid.coordinates()[j]] = depl_radial0(pos0[j],fluid.coordinates()[j],depl,pos_centre,L/2,dt)
vitesse[0] = liste_depl[k][0]/dt
vitesse[1] = liste_depl[k][1]/dt
v_part[0] = liste_depl[k][0]/dt
v_part[1] = liste_depl[k][1]/dt
print(v_part)
p.set_boundary_velocity("Grains", v_part)
fluid.set_open_boundary("Grains", velocity=vitesse)
k+=1
depl[0] += liste_depl[k][0]
depl[1] += liste_depl[k][1]
print(k)
if k==(n-1):
vitesse[0] = 0
vitesse[1] = 0
v_part[0] = 0
v_part[1] = 0
p.set_boundary_velocity("Grains", v_part)
fluid.set_open_boundary("Grains", velocity=vitesse)
k+=1
for j in range(N):
fluid.velocity_mesh()[j] = vitesse
fluid.implicit_euler(dt, newton_atol=5e-6, newton_rtol=1e-5, newton_max_it=150, reduced_gravity=0, stab_param=0.)
# Computation of the new velocities
forces = fluid.compute_node_force(dt)
vn = p.velocity() + forces*dt/p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
# Number of contact 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()))
# NLGS iterations
for i in range(nsub) :
tol = 1e-6 #numerical tolerance for particles intersection
p.iterate(dt/nsub, forces, tol)
# Localisation of the particles in the fluid
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())