betonlmgc.py

#!/usr/bin/env python

#variable to use lmgc or not
use_lmgc = 1


from marblesbag import fluid as fluid
from marblesbag import scontact2
if use_lmgc :
    from marblesbag import lmgc2Interface

import numpy as np
import os
import time
import shutil

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
rhop is the particles density'''
def genInitialPosition(filename, r, rout, rhop) :
    p = scontact2.ParticleProblem()
    #Loading of the mesh.msh file specifying physical boundaries name
    p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
    #Definition of the points where the grains are located
    x = np.arange(rout, -rout, 2 * -r)
    y = np.arange(-rout, -2*rout/3., 2 * r)
    x, y = np.meshgrid(x, y)
    R2 = x**2 + y**2
    #condition to be inside the outer boundary
    keep = R2 < (rout - r)**2
    x = x[keep]
    y = y[keep]
    #First layer of grains
    for i in range(x.shape[0]) :
        p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
    
    x = np.arange(rout, -rout, -r)
    y = np.arange(-2*rout/3.+r, -rout/2.,  r)
    x, y = np.meshgrid(x, y)
    R2 = x**2 + y**2
    keep = R2 < (rout - r)**2
    x = x[keep]
    y = y[keep]
    #Second layer of grains
    for i in range(x.shape[0]) :
        p.add_particle((x[i], y[i]), r/2., r**2/4. * np.pi * rhop);
        
    x = np.arange(rout, -rout, 2 * -r/3)
    y = np.arange(-rout/2.+r, -5*rout/12., 2 * r/3)
    x, y = np.meshgrid(x, y)
    R2 = x**2 + y**2
    keep = R2 < (rout - r/3)**2
    x = x[keep]
    y = y[keep]
    #Third layer of grains
    for i in range(x.shape[0]) :
        p.add_particle((x[i], y[i]), r/3, r**2 /9 * np.pi * rhop);
    p.write_vtk(filename,0,0)


t = 0
ii = 0

#physical parameters
g =  -9.81                              # gravity
rho = 1e3                               # fluid density
rhop = 1600                             # grains density
nu = 1e-4                               # kinematic viscosity
tEnd = 100                              # final time

#numerical parameters
h = 0.002                               # mesh size
dt = 5e-3                               # time step
epsilon = 1e-4                          # stabilization parametre

#Geometrical parametres
rout = 0.0254                           # outer boundary
r = 397e-6/2                            # grains radius
shutil.copy("mesh.msh", outputdir +"/mesh.msh")

#Object particles creation
genInitialPosition(outputdir, r, rout, rhop)
if use_lmgc :
    friction=0.1                        # friction coefficient
    lmgc2Interface.scontactToLmgc90(outputdir,0,friction)
    p = lmgc2Interface.LmgcInterface()
else :
    p = scontact2.ParticleProblem()
    p.read_vtk(outputdir,0)

outf = 5                                # number of iterations between output files
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Outer",0,lambda x : -4*x[:, 1]),("Outer",1,lambda x : 4*x[:, 0]),("Inner",0,0.),("Inner",1,0.),("PtFix",2,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())
fluid.export(outputdir,0,0)

tic = time.clock()
#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))