Commit b9f091ba authored by Nathan Coppin's avatar Nathan Coppin
Browse files

updating testcases part 2

parent 86e8b17c
Pipeline #8909 passed with stages
in 5 minutes and 40 seconds
......@@ -40,9 +40,9 @@ def genInitialPosition(filename, r, rout, rhop) :
rin -- the inner radius of the geometry
rhop -- the particles density
"""
p = scontact.ParticleProblem(2,True,True)
p = scontact.ParticleProblem(2,True)
#Loading of the mesh.msh file specifying physical boundaries name and material
p.load_msh_boundaries("mesh.msh", ["Outer","Top",],material="PVC")
p.load_msh_boundaries("mesh.msh", ["Outer"],material="PVC")
#Definition of the points where the grains are located
x = np.arange(rout-1.2*r, -rout+1.2*r, -2*r*1.05)
x, y = np.meshgrid(x, x)
......@@ -68,23 +68,21 @@ tEnd = 10 # final time
t = 0
#numerical parameters
dt = 1e-3 # time step
outf = 100 # number of iterations between output files
outf = 10 # number of iterations between output files
ii = 0
#geometry parameters
rout = 0.05 # outer radius
r = 1e-3 # grains radius
rout = 0.11 # outer radius
r = 2e-3 # grains radius
v = 0.1 # tangent velocity
#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, rout, rhop)
p = scontact.ParticleProblem(2,True,True)
p = scontact.ParticleProblem(2,True)
p.read_vtk(outputdir,0)
p.set_friction_coefficient(.3,"Glass","Glass") # between particles
p.set_friction_coefficient(.3,"Glass","PVC") # between particles and boundaries
p.set_tangent_boundary_velocity("Outer",-v) # setting the tangent velocity of the rotating mill
p.set_tangent_boundary_velocity("Top",-v) # setting the tangent velocity of the rotating mill
#
# FLUID PROBLEM
#
......@@ -92,10 +90,8 @@ fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Outer",0,lambda x : (x[:,1]/rout)*-v )
fluid.set_strong_boundary("Outer",1,lambda x : (-x[:,0]/rout)*-v)
fluid.set_wall_boundary("Top",0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.solution()[:,2] = g[1]*rho*(fluid.coordinates()[:,1] - rout/2)
fluid.export_vtk(outputdir, 0.,0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.write_vtk(outputdir, 0,0)
#
# COMPUTATION LOOP
#
......@@ -106,5 +102,5 @@ while t < tEnd :
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
# Migflow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of Marblesbag: 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
from migflow import scontact
from migflow import fluid
from migflow import lmgc90Interface
from migflow import time_integration
from pylmgc90 import pre
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, rhop) :
p = scontact.ParticleProblem(3,True)
#Loading of the mesh.msh file specifying physical boundaries name and material
p.load_msh_boundaries("mesh.msh", ["Outer","Front","Rear"],material="PVC")
# Positioning the particles on a regular grid
for y in np.arange(.04-r*1.11, -0.04+r*1.11, -2*r*1.03):
for z in np.arange(.0+r*1.11, .01-r*1.11, 2*r*1.03):
for x in np.arange(.1-r*1.11, -.1+r*1.11, -2*r*1.03):
if x**2 + y**2 <= (rout-2*r)**2:
r1 = np.random.normal(0.9*r,0.05*r)
p.add_particle((x, y, z), r1, 4./3.* r1**3 * np.pi * rhop,"Glass");
p.write_vtk(filename,0,0)
t = 0
ii = 0
#physical parameters
g = np.array([0,-9.81,0]) # gravity
rhop = 2500 # grains density
rho = 1000
nu = 1e-6
tEnd = 10 # final time
#numerical parameters
dt = 1e-3 # time step
#geometry parameters
rout = 0.05 # outer radius
r = 1.6e-3 # grains radius
Fr = 4e-1 #Froude number
v = (Fr*rout*9.81)**0.5 #Corresponding tangent velocity
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#Object particles creation
genInitialPosition(outputdir, r, rout, rhop)
p = scontact.ParticleProblem(3,True)
p.read_vtk(outputdir,0)
outf = 10 # number of iterations between output files
#Taking friction into account
p.set_friction_coefficient(.2,"Glass","Glass")#Between particles
p.set_friction_coefficient(.2,"Glass","PVC")#Between particles and boundaries
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Outer",0,lambda x : (x[:,1]/rout)*-v )
fluid.set_strong_boundary("Outer",1,lambda x : (-x[:,0]/rout)*-v)
fluid.set_strong_boundary("Front",0,lambda x : (x[:,1]/rout)*-v*((x[:,1]**2+x[:,0]**2)**0.5)/rout)
fluid.set_strong_boundary("Front",1,lambda x : (-x[:,0]/rout)*-v*((x[:,1]**2+x[:,0]**2)**0.5)/rout)
fluid.set_strong_boundary("Rear",0,lambda x : (x[:,1]/rout)*-v*((x[:,1]**2+x[:,0]**2)**0.5)/rout)
fluid.set_strong_boundary("Rear",1,lambda x : (-x[:,0]/rout)*-v*((x[:,1]**2+x[:,0]**2)**0.5)/rout)
# 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())
fluid.write_vtk(outputdir, 0,0)
while t < tEnd :
time_integration.iterate(fluid,p,dt,min_nsub=5,contact_tol=1e-6,external_particles_forces=g*p.mass())
t += dt
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
# 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
# Bidimensional particles sedimentation in fluid
from migflow import time_integration
from migflow import fluid
from migflow import scontact
import numpy as np
import os
import time
import shutil
import random
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameter
g = np.array([0,-9.81]) # gravity
r = .1 # particles radius
rhop = 1500 # particles density
rho = 800 # fluid density
nu = 0.9e-6 # kinematic viscosity
mf = 10
# Numerical parameters
outf = 5 # number of iterations between output files
dt = 1e-3 # time step
tEnd = 20 # final time
# Geometrical parameters
lx = 3.5
#
# PARTICLE PROBLEM
#
p = scontact.ParticleProblem(2,True)
p.load_msh_boundaries("mesh.msh", ["Bottom", "Lateral", "Top"])
p.add_particle((0, 0), r, mf,"Name", forced=1)
p.add_particle((0, .25), r, mf,"Name",forced=1)
p.add_particle((0, .5), r, mf,"Name",forced=1)
p.add_particle((0, .75), r, mf,"Name",forced=1)
p.add_particle((.25, .75), r, mf,"Name",forced=1)
p.add_particle((.5, .75), r, mf,"Name",forced=1)
p.add_particle((.25, .25), r, mf,"Name",forced=1)
p.position()[:,0] -= 0.85
p.add_particle((0, 0), r, mf,"Name",forced=1)
p.add_particle((0, .25), r, mf,"Name",forced=1)
p.add_particle((0, .5), r, mf,"Name",forced=1)
p.add_particle((0, .75), r, mf,"Name",forced=1)
p.add_particle((.25, 0), r, mf,"Name",forced=1)
p.add_particle((.5, 0), r, mf,"Name",forced=1)
p.position()[:,0] -= 0.85
p.add_particle((0, 0), r, mf,"Name",forced=1)
p.add_particle((0, .25), r, mf,"Name",forced=1)
p.add_particle((0, .5), r, mf,"Name",forced=1)
p.add_particle((0, .75), r, mf,"Name",forced=1)
p.add_particle((.25, .75), r, mf,"Name",forced=1)
p.add_particle((0.5, 0), r, mf,"Name",forced=1)
p.add_particle((0.5, .25), r, mf,"Name",forced=1)
p.add_particle((0.5, .5), r, mf,"Name",forced=1)
p.add_particle((0.5, .75), r, mf,"Name",forced=1)
p.add_particle((.25, 0), r, mf,"Name",forced=1)
p.position()[:,0] -= 0.85
p.add_particle((0, 0), r, mf,"Name",forced=1)
p.add_particle((0, .25), r, mf,"Name",forced=1)
p.add_particle((0, .5), r, mf,"Name",forced=1)
p.add_particle((0, .75), r, mf,"Name",forced=1)
p.add_particle((.25, .25), r, mf,"Name",forced=1)
p.add_particle((0.5, 0.5), r, mf,"Name",forced=1)
p.add_particle((0.75, .25), r, mf,"Name",forced=1)
p.add_particle((1, .5), r, mf,"Name",forced=1)
p.add_particle((1, .75), r, mf,"Name",forced=1)
p.add_particle((1, 0.25), r, mf,"Name",forced=1)
p.add_particle((1, 0), r, mf,"Name",forced=1)
p.position()[:,1] -= 1.1
p.add_particle((0, 0), r, mf,"Name",forced=1)
p.add_particle((0, .25), r, mf,"Name",forced=1)
p.add_particle((0, .5), r, mf,"Name",forced=1)
p.add_particle((0, .75), r, mf,"Name",forced=1)
p.add_particle((.25, .5), r, mf,"Name",forced=1)
p.add_particle((0.5, 0.25), r, mf,"Name",forced=1)
p.add_particle((0.75, .5), r, mf,"Name",forced=1)
p.add_particle((1, .5), r, mf,"Name",forced=1)
p.add_particle((1, .75), r, mf,"Name",forced=1)
p.add_particle((1, 0.25), r, mf,"Name",forced=1)
p.add_particle((1, 0), r, mf,"Name",forced=1)
p.position()[p.position()[:,1]>=0,0] -= 1.35
p.add_particle((0, 0), r, mf,"Name",forced=1)
p.add_particle((0, .25), r, mf,"Name",forced=1)
p.add_particle((0, .5), r, mf,"Name",forced=1)
p.add_particle((0, .75), r, mf,"Name",forced=1)
p.position()[p.position()[:,1]>=0,0] -= 0.35
p.add_particle((0, 0), r, mf,"Name",forced=1)
p.add_particle((0, .25), r, mf,"Name",forced=1)
p.add_particle((0, .5), r, mf,"Name",forced=1)
p.add_particle((0, .75), r, mf,"Name",forced=1)
p.add_particle((0.25, 0), r, mf,"Name",forced=1)
p.add_particle((0.5, 0), r, mf,"Name",forced=1)
p.add_particle((0.75, 0), r, mf,"Name",forced=1)
p.add_particle((0.75, .25), r, mf,"Name",forced=1)
p.add_particle((0.75, 0.5), r, mf,"Name",forced=1)
p.add_particle((0.5, .5), r, mf,"Name",forced=1)
p.add_particle((0.25, .75), r, mf,"Name",forced=1)
p.add_particle((0.5, .75), r, mf,"Name",forced=1)
p.add_particle((0.75, .75), r, mf,"Name",forced=1)
p.position()[p.position()[:,1]>=0,0] -= 0.35
p.position()[:,0] += 0.8625
r = 0.023
x = np.arange(-3+r, 3-r, 2.5*r)
y = np.arange(4-r, 1 +r, -2.5*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,"Sand")
p.set_friction_coefficient(0.4,"Sand","Name")
p.set_friction_coefficient(0.2,"Sand","Sand")
p.write_vtk(outputdir,0,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("Top",pressure=0)
fluid.set_wall_boundary("Lateral", velocity = [0,0])
fluid.set_wall_boundary("Bottom")
# 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(),np.zeros_like(p.velocity()))
fluid.export_vtk(outputdir,0,0)
tic = time.time()
while t < tEnd :
time_integration.iterate(fluid,p,dt,external_particles_forces=g*p.mass())
t += dt
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 4000
p.write_vtk(outputdir, ioutput, t)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......@@ -84,7 +84,7 @@ ii = 0
t = ii*dt
#set initial_condition
fluid.export_vtk(outputdir,t,ii)
fluid.write_vtk(outputdir,ii,t)
tic = time.time()
while t < tEnd:
......@@ -94,7 +94,7 @@ while t < tEnd:
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......
# MigFlow - Copyright (C) <2010-2020>
# <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
# Deposit of 3D particles in a box full of fluid.
from migflow import fluid
from migflow import scontact
from migflow import time_integration
import numpy as np
import os
import time
import shutil
def genInitialPosition(p, filename, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure.
For this testcase, radius are contained in a file radius.csv created with a special distribution
Keyword arguments:
filename -- name of the output file
rhop -- particles density
"""
# Reading the file containing the radii
myRadius = np.genfromtxt('radius.csv', delimiter=',')*2
r2 = np.amax(myRadius)*8e-6
i = 0
# Positioning the particles on a regular grid
for y in np.arange(0, -.05, -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]*8e-6, 4./3.* (myRadius[i%500000]*8e-6)**3 * np.pi * rhop);
i += 1
# Translation of the particles to the top of the box
state = p.state()
state.x[:, 1] += .07/2
p.set_state(state)
print("Number of particles=%d" % len(p.state().x[:, 1]))
# p.write_vtk(filename, 0, 0)
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
# g = np.array([0,-9.81,0]) # gravity
g = np.array([0,0,0])
rho = 1000 # fluid density
rhop = 2640 # particles density
nu = 1e-6 # kinematic viscosity
L = 1
# Numerical parameters
outf = 1 # number of iterations between output files
dt = 0.1 # time step
tEnd = 2 # final time
#
# PARTICLE PROBLEM
#
# Particles structure creation and setting particles properties
r = 1e-2
p = scontact.ParticleProblem(3)
p.load_msh_boundaries("mesh.msh", ["Bottom", "Top", "X", "Z"])
p.add_particle((r,0.5,0.5), r, 4/3*r**3*np.pi*rhop)
p.add_particle((L-2*r,0.5,0.5), r, 4/3*r**3*np.pi*rhop)
state = p.state()
state.v[1, 0] += 0.5
p.set_state(state)
# genInitialPosition(p,outputdir, rhop)
p.write_vtk(outputdir,0,0)
# Initial time and iteration
t = 0
ii = 0
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho])
#Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("X")
fluid.set_wall_boundary("Z")
#Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces(),init=True)
tic = time.time()
fluid.export_vtk(outputdir,0,0)
G = np.zeros((p.n_particles(),p.dim()))
G[:,1] = p.mass()[:,0]*g[1]
#
# COMPUTATION LOOP
#
while t < tEnd :
time_integration.iterate(fluid, p, dt,min_nsub=20, contact_tol = 1e-8,external_particles_forces=G)
state = p.state()
state.x[:,0] = np.fmod(state.x[:,0],L)
ind = np.where(state.x[:,0] <= 0)
state.x[ind,0] += L
p.set_state(state)
t += dt
# Output files writing
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))
This diff is collapsed.
......@@ -62,7 +62,7 @@ fluid.set_wall_boundary("Top",velocity=[0,0])
fluid.set_open_boundary("RightUp",pressure=0)
fluid.set_open_boundary("RightDown",pressure=0)
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
tic = time.time()
#
# COMPUTATION LOOP
......@@ -78,7 +78,7 @@ while t < tEnd :
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
s = fluid.solution()
......@@ -86,4 +86,4 @@ x = fluid.coordinates()
vel = (s[:,0]-1/(20*nu*rho)*x[:,1]*(1-x[:, 1]))**2
vS = np.sum((1/(20*nu*rho)*x[:,1]*(1-x[:, 1]))**2)
print(f"Error criteria : {(vel.sum())**.5} \t < {(vS**0.5)/50}")
print(f"Succeed ? {(vel.sum())**.5 < (vS**0.5)/50}")
\ No newline at end of file
print(f"Succeed ? {(vel.sum())**.5 < (vS**0.5)/50}")
#!/usr/bin/env python
use_lmgc90 = False #True
from migflow import fluid
from migflow import scontact
if use_lmgc90 :
from migflow import lmgc90Interface
from migflow import time_integration
import numpy as np
import os
......@@ -33,7 +29,7 @@ t = 0
ii = 0
#physical parameters
g = 0. #-9.81
g = np.array([0,0])
rho = 1200
rhop = 1800 #2400
nu = 1e-4
......@@ -50,49 +46,30 @@ y0part = 0
lypart = 150*r
lxpart = 50*r
genInitialPosition(outputdir, r, x0part, y0part, lxpart, lypart, rhop)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
if use_lmgc90 :
# friction coefficient
friction = 0.1
lmgc90Interface.scontactTolmgc90(outputdir,2,0,friction)
p = lmgc90Interface.ParticleProblem(2)
else :
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
deltamass = p.volume().sum()*(rhop-rho)
deltapinit = deltamass*-g/lxpart
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 20 #000
outf = 10 #000
fluid = fluid.FluidProblem(2,g,nu*rho,rho)
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Box",0,0)
fluid.set_strong_boundary("Box",1,0)
fluid.set_open_boundary("Left",velocity=[0.001, 0])
fluid.set_open_boundary("Left",velocity=[0.0005, 0])