depot.py

#!/usr/bin/env python
from marblesbag import fluid3 as fluid
from marblesbag import scontact3

import numpy as np
import os
import time
import shutil

#grains creation and placement + physical boundaries definition
#The radii of the grains are defined in radius.csv and are computed to fit the distribution size given in the article of Andre et al. (2011) "Real-time analysis of the growth of granular media by an ultrasonic method: Application to the sedimentation of glass balls in water"
def genInitialPosition(filename, rhop) :
    p = scontact3.ParticleProblem()
    #Loading of the mesh.msh file specifying physical boundaries name
    p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
    #Reading the file containing the radius
    myRadius = np.genfromtxt('radius.csv', delimiter=',')
    r2=np.amax(myRadius)*1e-6
    i = 0
    #Positioning the grains 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
    p.position()[:, 1] += .15/2 
    print("Number of grains=%d" % len(p.position()[:, 1]))
    p.write_vtk(filename, 0, 0)


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

t = 0
ii = 0

p = scontact3.ParticleProblem()

#physical parameters
g =  -9.81                                      # gravity
rho = 1000                                      # fluid density
rhop = 2640                                     # grains density
nu = 1e-6                                       # kinematic viscosity
tEnd = 100                                      # final time

#numerical parameters
lcmin = 0.001                                   # mesh size
dt = 1e-3                                       # time step
alpha = 2.5e-4                                  # stabilization coefficient
epsilon = alpha*lcmin**2 /nu                    # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")

#Object particles creation
genInitialPosition(outputdir, rhop)
p = scontact3.ParticleProblem()
p.read_vtk(outputdir,0)

outf = 100                                      # number of iterations between output files

#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal-width velocity; field 1 is vertical velocity; field 2 is the horizontal-depth velocity; field 3 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

tic = time.clock()
fluid.export_vtk(outputdir,0,0)
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
    #Fluid solver
    fluid.implicit_euler(dt)
    forces = fluid.compute_node_force(dt)
    #Computation of the new velocities
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    #number of 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
    for i in range(nsub) :
        tol = 1e-6
        p.iterate(dt/nsub, forces, tol)
    #Localisation of the grains in the fluid
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    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.clock() - tic))