drop.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
import random

def genInitialPosition(filename, r, rout, rhop, compacity) :
    p = scontact3.ParticleProblem()
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom","X","Z"])
    
    N = compacity*rout**3/(r**3)
    e = 2*rout/(6*N/np.pi)**(1./3.)
    
    for x in np.arange(rout, -rout, -e):
        for y in np.arange(rout, -rout, -e):
            for z in np.arange(rout, -rout, -e):
                R2 = x**2 + y**2 + z**2
                if R2<(rout-r)**2:
                    p.add_particle((x+random.uniform(-e/2+r,e/2-r), y+random.uniform(-e/2+r,e/2-r), z+random.uniform(-e/2+r,e/2-r)), r, 4./3.* r**3 * np.pi * rhop);
    
    p.position()[:, 1] += 0.52
    print(len(p.position()[:, 1] ))
    p.write(filename)


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

t = 0
ii = 0

r=154e-6
R = 3.3e-3
compacity = .20
drho = 35.4
p = scontact3.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))

#physical parameters
g =  -9.81
rhop = 2450
nu = 1.17/1030.
rho = 1030#rhop-drho/compacity
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 1000

#numerical parameters
lcmin = 0.0004 # approx r*100 but should match the mesh size
dt = 5e-5*1000
alpha = 1e-3
epsilon = alpha*lcmin**2 /nu
print('epsilon',epsilon)

shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
p.write(outputdir+"/part-00000")
mu = nu*rho

genInitialPosition(outputdir + "/part-00000", r, R, rhop, compacity)

p = scontact3.ParticleProblem(outputdir+"/part-00000")

print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 10
outf1 = 100000

#strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Box",1,0.),("Box",0,0.),("Box",2,0.),("Top",3,0.)]
#strong_boundaries = [("Top",0,.0),("Top",1,0.),("Top",2,0.),("Bottom",0,.0),("Bottom",1,0.),("Bottom",2,0.),("X",0,.0),("X",1,0.),("X",2,0.),("Z",0,.0),("Z",1,0.),("Z",2,0.),("Top",3,0.)]
strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("X",0,.0),("Z",2,0.),("Top",3,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries, 1)

fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())


ii = 0
jj = 0
tEnd1 = 100
#p.velocity()[:,1] = 2./9. *(compacity)* R**2*g*(rhop-rho)/mu
print(p.velocity()[1,1])
print(1-compacity)


for jj in range(4):
    dt1 = dt
    t = 0
    if jj!=0:
        fluid.adapt_mesh(1e-3,1e-4,300000, epsilon, 1)
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    
    while t < tEnd1 : 
        niter = fluid.implicit_euler(dt1)
        forces = fluid.compute_node_force(dt1)
        vn = p.velocity() + forces * dt1 / p.mass()
        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
        if niter < 5 and dt1 < 20:
            dt1 *= 3
        elif niter > 5 and dt1 > dt:
            dt1 /= 3
        print(dt1)
        t += dt1
#        if ii %outf == 0 :
#                ioutput = int(ii/outf) + 1
#                p.write_vtk(outputdir, ioutput, t)
#                fluid.export_vtk(outputdir, t, ioutput)
        ii += 1

ii = 0
t = 0
tic = time.clock()
fluid.export_vtk(outputdir,0,0)
forces = g*p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd : 
    if (ii%5==0 and ii != 0):
       fluid.adapt_mesh(1e-3,1e-4,300000, epsilon, 1)
       fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    fluid.implicit_euler(dt)
    forces = fluid.compute_node_force(dt)
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    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()))
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
    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))