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

depot 2d/3d

parent 8cd8c59f
Pipeline #4501 failed with stage
in 13 seconds
Hide whitespace changes
Inline Side-by-side
testcases/depot-2d/dep.py → testcases/depot-2d/depot.py
View file @ acee8278
......@@ -21,11 +21,10 @@
#!/usr/bin/env python
'''
Mixing of grains in a rotating drum.
This example shows how to set boundary conditions as a function of some parametres.
A boolean parametre give the possibility to use lmgc90 to solve contacts or not.
'''
# TESTCASE DESCRIPTION
# Bidimensional particles sedimentation in fluid
from migflow import fluid
from migflow import scontact
......@@ -35,130 +34,117 @@ import time
import shutil
import random
'''
genInitialPosition is a function that sets all the grains centre positions and create the grains objects to add in the computing structure
Args: -filename : name of the output file
-r : max radius of the grains
-H : domain height
-ly : grains area height
-lx : grains area width
-rhop : grains density
'''
def genInitialPosition(filename, r, H, ly, lx, rhop) :
#Grains structure builder
p = scontact.ParticleProblem(2)
#Load mesh.msh file specifying physical boundaries names
"""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"])
#Definition of the points where the grains are located
x = np.arange(-lx/2+r, lx/2-r, 2*r)
y = np.arange(H/2-r, H/2-ly+r, -2*r)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
#Add a grain at each centre position
#Definition of the points where the particles are located
x = np.arange(-lx/2+r, lx/2-r, 2*r)
y = np.arange(H/2-r, H/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
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
#physical parameters
g = -9.81 # gravity
r = 1e-3 # grains radius
rhop = 1500 # grains density
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
#numerical parameters
outf = 5 # number of iterations between output files
dt = 1e-2 # time step
tEnd = 100 # final time
#geometrical parameters
ly = 5e-2 # grains area height
lx = 4e-1 # grains area widht
H = 1.2 # domain height
#Initialise grains
# Physical parameters
g = -9.81 # gravity
r = 1e-3 # particles radius
rhop = 1500 # particles density
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
# Numerical parameters
outf = 5 # number of iterations between output files
dt = 1e-2 # time step
tEnd = 100 # final time
# Geometrical parameters
ly = 5e-2 # particles area height
lx = 4e-1 # particles area widht
H = 1.2 # domain height
#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, H, ly, lx, rhop)
p = scontact.ParticleProblem(2)
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
# Initial time and iteration
t = 0
ii = 0
'''
fluid_problem is the function builder of the fluid class
Args: -g : gravity
-viscosity : list of fluid phases viscosities
-density : list of fluid phases densities
'''
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
#Set the mesh geometry for the fluid computation
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
'''
set_strong_boundary is the function setting the strong boundaries (=constraints) for the fluid problem
Args: -Tag : physical tag (set in the mesh.geo file) of the boundary on which you want to add a constraint
-Field : O: x-velocity; 1: y-velocity; 2: pressure
-Value : value assigned to the specified field on the specified boundary
'''
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Lateral",0,0)
'''
set_weak_boundary is the function setting the weak boundaries (=normal fluxes) for the fluid problem
Args: -Tag : physical tag (set in the mesh.geo file) of the boundary on which you want to add a constraint
-Bnd Type : boundaries can be of type: wall, inflow, outflow
'''
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Wall")
#Set location of the grains in the mesh and compute the porosity in each computation cell
# 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())
fluid.export_vtk(outputdir,0,0)
tic = time.time()
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#
# COMPUTATION LOOP
#
while t < tEnd :
#Fluid solver
# Fluid solver
fluid.implicit_euler(dt)
#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()))))
# 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()))
#Contact solver
# NLGS iterations
for i in range(nsub) :
tol = 1e-6 #numerical tolerance for grains intersection
tol = 1e-6 #numerical tolerance for particles intersection
p.iterate(dt/nsub, forces, tol)
#Localisation of the grains in the fluid
# Localisation of the particles in the fluid
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
t += dt
t += dt
#Output files writting
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
testcases/depot-3d/depot.py
View file @ acee8278
......@@ -20,80 +20,81 @@
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
from migflow import fluid3 as fluid
from migflow import scontact3
# TESTCASE DESCRIPTION
# Deposit of 3D particles in a box full of fluid.
from migflow import fluid
from migflow import scontact
import numpy as np
import os
import time
import shutil
'''
Deposit of 3D grains in a box full of fluid.
'''
'''
genInitialPosition is a function that sets all the grains centre positions and create the grains objects to add in the computing structure.
For this testcase, radius are contained in a file radius.csv created with a special distribution
Args: -filename : name of the output file
-rhop : grains density
'''
def genInitialPosition(filename, rhop) :
#Grains structure builder
p = scontact3.ParticleProblem()
#Load the mesh.msh file specifying physical boundaries names
"""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
"""
# Particles structure builder
p = scontact.ParticleProblem(3)
# Load the mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
#Reading the file containing the radii
# Reading the file containing the radii
myRadius = np.genfromtxt('radius.csv', delimiter=',')
r2 = np.amax(myRadius)*1e-6
i = 0
#Positioning the grains on a regular grid
r2 = np.amax(myRadius)*1e-6
i = 0
# Positioning the particles 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
# Translation of the particles to the top of the box
p.position()[:, 1] += .15/2
print("Number of grains=%d" % len(p.position()[:, 1]))
print("Number of particles=%d" % len(p.position()[:, 1]))
p.write_vtk(filename, 0, 0)
#Define output directory
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#physical parameters
g = -9.81 # gravity
rho = 1000 # fluid density
rhop = 2640 # grains density
nu = 1e-6 # kinematic viscosity
# Physical parameters
g = -9.81 # gravity
rho = 1000 # fluid density
rhop = 2640 # particles density
nu = 1e-6 # kinematic viscosity
#numerical parameters
outf = 100 # number of iterations between output files
dt = 1e-3 # time step
tEnd = 100 # final time
# Numerical parameters
outf = 100 # number of iterations between output files
dt = 1e-3 # time step
tEnd = 100 # final time
#Grains structure creation and setting grains properties
#
# PARTICLE PROBLEM
#
# Particles structure creation and setting particles properties
genInitialPosition(outputdir, rhop)
p = scontact3.ParticleProblem()
p = scontact.ParticleProblem(3)
p.read_vtk(outputdir,0)
'''
fluid_problem is the function builder of the fluid class
Args: -g : gravity
-viscosity : list of fluid phases viscosities
-density : list of fluid phases densities
'''
fluid = fluid.fluid_problem(g,[nu*rho],[rho])
# 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")
'''
set_strong_boundary is the function setting the strong boundaries (=constraints) for the fluid problem
Args: -Tag : physical tag (set in the mesh.geo file) of the boundary on which you want to add a constraint
-Field : O: x-velocity; 1: y-velocity; 2: pressure
-Value : value assigned to the specified field on the specified boundary
'''
fluid.set_strong_boundary("Top",0,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Top",2,0)
......@@ -101,40 +102,41 @@ fluid.set_strong_boundary("Top",3,0)
fluid.set_strong_boundary("Box",0,0)
fluid.set_strong_boundary("Box",1,0)
fluid.set_strong_boundary("Box",2,0)
#Set location of the grains in the mesh and compute the porosity in each computation cell
#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())
fluid.export_vtk(outputdir,0,0)
tic = time.clock()
#Computation loop
tic = time.clock()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#Computation loop
#
# COMPUTATION LOOP
#
while t < tEnd :
#Fluid solver
# Fluid solver
fluid.implicit_euler(dt)
#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
# Adaptation of the mesh.
if (ii%15==0 and ii != 0):
fluid.adapt_mesh(1e-2,3e-3,10000)
#Computation of the new velocities
# 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()))))
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()))
#Contact solver
# NLGS iterations
for i in range(nsub) :
p.iterate(dt/nsub, forces)
t += dt
t += dt
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#Output files writing
# 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))
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
\ No newline at end of file
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