shaker.py

# 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
# Sort particles with respect to their mass by jigging the granular matrix
from migflow import fluid
from migflow import scontact
from migflow import time_integration

import numpy as np
import os
import time
import shutil
import random

def genInitialPosition(filename, r, ly, 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 -- radius of the particles
    ly - particles area height
    rhop -- particles density        
    """
    p = scontact.ParticleProblem(2)
    # Loading of the mesh.msh file specifying physical boundaries name
    p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom","BottomOut","TopOut"])
    # Definition of the points where the grains are located
    x = np.arange(-0.5+r, 0.5-r+1e-10, 2*r)
    y = np.arange(-ly/2+r, 0.5*ly/2-r, 2*r)
    x, y = np.meshgrid(x, y)
    x = x.flat
    y = y.flat
    for i in range(len(x)) :
        rhop = random.choice([1100,1350,1600])
        # Addition of a particle object at each point
        p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
    p.write_vtk(filename,0,0)

outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

# physical parameters
g =  -9.81                                  # gravity
rho = 1000                                  # fluid density
rhop = 1500                                 # grains density
nu = 1e-6                                   # kinematic viscosity
r = 5e-3                                    # grains radius
ly = 1                                      # box height

# numerical parameters
dt = 5e-3                                   # time step
tEnd = 100                                  # final time
outf = 10                                   # number of iterations between output files

#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, ly, rhop)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)

print("r = %g, m = %g" %  (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
print("number of grains = %d"%len(p.r()))

# 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")
# periodic boundary condition function
def outerBndV(x) :
    return 0.15*max(np.sin(-t*np.pi*2./5-4*np.pi/5),0)
fluid.set_wall_boundary("Lateral")
fluid.set_open_boundary("Bottom",velocity=[0, outerBndV])
fluid.set_open_boundary("BottomOut",velocity=[0, outerBndV])
fluid.set_open_boundary("TopOut",pressure=0)
fluid.set_open_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(), init=True)
fluid.export_vtk(outputdir,0,0)

tic = time.time()
G = np.zeros_like(p.velocity())
G[:,1] = p.mass()[:,0]*g
#
# COMPUTATION LOOP
#
while t < tEnd :
    time_integration.iterate(fluid, p, dt, min_nsub=25, external_particles_forces=G)
    t += dt
    # 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
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))