Commit 67fb6938 authored by Michel Henry's avatar Michel Henry
Browse files

Granular Collapse

parent 8d2d1d82
Pipeline #8782 passed with stages
in 5 minutes and 18 seconds
......@@ -74,6 +74,7 @@ def showFluidVelocity(fluid, renderView):
fluidVelocity.Function = 'velocity/porosity'
RenameSource('Fluid Velocity', fluidVelocity)
fluidVelocityDisplay = Show(fluidVelocity, renderView)
ColorBy(fluidVelocityDisplay, ('POINTS', 'FluidVelocity', 'Magnitude'))
fluidVelocityDisplay.SetScalarBarVisibility(renderView, True)
Hide(fluidVelocity)
return fluidVelocity
......
......@@ -433,10 +433,13 @@ class ParticleProblem :
x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
tagnames = VTK.string_array_decode(fdata["TagNames"])
tagmap = {}
print(fdata)
print(tagnames)
for j,n in enumerate(tagnames) :
tagmap[j] = self._lib.particleProblemGetTagId(self._p,n)
offsets = np.hstack([[0],cells["offsets"]])
el = cells["connectivity"]
print(tagmap)
tags = np.vectorize(tagmap.get)(cdata["Tag"])
materials = cdata["Material"]
for idx in np.where(cells["types"] == 1)[0] :
......
from migflow import time_integration
from migflow import scontact
import numpy as np
import os
import subprocess
import time
import shutil
import random
SEED = np.random.seed(0)
def genInitialPosition(filename, r, H, 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
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("meshDepot.msh", ["Top", "Lateral", "Bottom"])
#Definition of the points where the particles are located
x = np.arange(r-1e-8, lx-r+1e-8, 2*r)
y = np.arange(r-1e-8, H-r-1e-8, 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)) :
R = r*(1-np.random.random()/2)
Xi = x[i] + np.random.random()/10000
p.add_particle((Xi, y[i]), R, R**2 * np.pi * rhop)
p.write_vtk(filename,0,0)
dir_path = os.path.dirname(os.path.realpath(__file__))
os.chdir(dir_path)
outputdir = "outputDepot"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
subprocess.call(["gmsh", "-2", "meshDepot.geo","-clscale","1"])
# Physical parameters
g = np.array([0,-9.81])
rhop = 9000 # particle density
r = .2e-2 # particle radius
# Geometrical parameters
L = 0.1
H = 0.1
# Numerical parameters
tEnd = 0.5
dt = 0.0025 # time step
outf = 5 # number of iterations between output files
genInitialPosition(outputdir,r,H,L,rhop)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
G = np.zeros((p.n_particles(),p.dim()))
G[:,1] = g[1]*p.mass()[:,0]
t = 0
ii = 0
tic = time.time()
#
# COMPUTATION LOOP
#
while t < tEnd :
# Time integration
time_integration.iterate(None,p,dt,external_particles_forces=G)
#Output files writting
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
\ No newline at end of file
# 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:
# Collapsing of a grain column in fluid.
from migflow import fluid
from migflow import scontact
from migflow import time_integration
import numpy as np
import os
import subprocess
import time
import shutil
import random
outputdir = "outputCollapse"
outputdirDepot = "outputDepot"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","1"])
# Physical parameters
g = np.array([0,-9.81]) # gravity
rhop = 9000 # grains density
rho = 1000 # fluid density
nu = 1e-6 # kinematic fluid viscosity
# Geometrical parameters
H = 0.15 # domain height
L = 0.3
# Numerical parameters
outf = 10 # number of iterations between output files
dt = 1e-3 # time step
tEnd = 10 # final time
#
# PARTICLES PROBLEM
#
# Particles container for the current problem
p = scontact.ParticleProblem(2)
p.load_msh_boundaries("mesh.msh", ["Top", "Lateral", "Bottom"])
pdepot = scontact.ParticleProblem(2)
pdepot.read_vtk(outputdirDepot,40)
p.add_particles(pdepot.position(),pdepot.r(),
pdepot.mass(),pdepot.velocity())
p.write_vtk(outputdir, 0, 0)
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho], [rho])
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom", velocity=[0,0])
fluid.set_wall_boundary("Lateral", velocity=[0,0])
fluid.set_wall_boundary("Top", velocity=[0,0], pressure=0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Lateral",0,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
fluid.solution()[:,2] = (H-fluid.coordinates()[:,1])*(-g[1])*rho
fluid.read_vtk(outputdir,0)
# Initial time and iteration.
t = 0
ii = 0
tic = time.time()
G = np.zeros((p.n_particles(),p.dim()))
G[:,1] = g[1]*p.mass()[:,0]
#
# COMPUTATION LOOP
#
while t < tEnd :
time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=G)
t += dt
# Output files writing.
if ii%outf == 0 :
ioutput = int(ii/outf) + 1
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))
L = 0.3;
H = 0.15;
y = 0;
lc = 0.01;
Point(1) = {0, 0, 0, lc};
Point(2) = {L, 0, 0, lc};
Point(3) = {L, H, 0, lc};
Point(4) = {0, H, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Point(5) = {0,0.5,0};
Line Loop(1) = {1:4};
Plane Surface(1) = {1};
Physical Line("Bottom") = {1};
Physical Line("Lateral") = {2,4};
Physical Line("Top") = {3};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
Mesh.CharacteristicLengthFromPoints = 1;
Mesh.MshFileVersion = 4;
L = 0.1;
H = 0.1;
y = 0;
lc = 0.01;
Point(1) = {0, 0, 0,lc};
Point(2) = {L, 0, 0,lc};
Point(3) = {L, H, 0,lc};
Point(4) = {0, H, 0,lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Point(5) = {0,0.5,0};
Line Loop(1) = {1:4};
Plane Surface(1) = {1};
Physical Line("Bottom") = {1};
Physical Line("Lateral") = {2,4};
Physical Line("Top") = {3};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
Mesh.CharacteristicLengthFromPoints = 1;
Mesh.MshFileVersion = 4;
......@@ -69,7 +69,6 @@ p.write_vtk(outputdir,0,0)
#Initial time and iteration
t = 0
ii = 0
ii = 0
tic = time.time()
......
Markdown is supported
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