depot-fluid.py 2.56 KB
Newer Older
dubois's avatar
maj  
dubois committed
1
#!/usr/bin/env python
2
3
from marblesbag import fluid as fluid
from marblesbag import scontact2
Matthieu Constant's avatar
Matthieu Constant committed
4
from marblesbag import lmgc2Interface
5

dubois's avatar
maj  
dubois committed
6
7
8
import numpy as np
import os
import time
9
import shutil
dubois's avatar
maj  
dubois committed
10
11
import random

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
12
def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
dubois's avatar
maj  
dubois committed
13
    p = scontact2.ParticleProblem()
14
    p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
15
16
    x = np.arange(ox+r, ox+lx, 2*r)
    y = np.arange(oy+r, oy+ly, 2*r)
dubois's avatar
maj  
dubois committed
17
18
    for i in range(x.shape[0]) :
       for j in range(y.shape[0]) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
19
          rr=r*(0.95+random.random()*0.05)
20
          p.add_particle((x[i], y[j]), rr, rr**2 * np.pi * rhop);
dubois's avatar
maj  
dubois committed
21
22
23
    p.write(filename)


24
outputdir = "output"
dubois's avatar
maj  
dubois committed
25
26
27
28
29
30
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

31
r=1e-4
dubois's avatar
maj  
dubois committed
32
33
34
p = scontact2.ParticleProblem()

#physical parameters
35
g =  -9.81
dubois's avatar
maj  
dubois committed
36
rho = 1200
37
rhop = 1800 #2400
Matthieu Constant's avatar
Matthieu Constant committed
38
nu = 1e-6
39
tEnd = 7
dubois's avatar
maj  
dubois committed
40
41

#numerical parameters
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
42
#h = 2e-4 # approx r*100 but should match the mesh size
dubois's avatar
maj  
dubois committed
43

44
dt = 1e-3 #0.5*r/V
Matthieu Constant's avatar
Matthieu Constant committed
45
epsilon = 1e-4 #2e-2*c*h # ?? not sure ??
dubois's avatar
maj  
dubois committed
46

Matthieu Constant's avatar
Matthieu Constant committed
47
48
print("eps = %g" % ( epsilon))
use_lmgc = True
49
50
y0part = 40e-3
lypart = 20e-3
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
51
52
lxpart = 0.015
genInitialPosition(outputdir + "/part-00000", r, 0., y0part, lxpart, lypart, rhop)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
53
54
55
56
57
if use_lmgc :
    friction=0.1
    lmgc2Interface.scontactToLmgc90(outputdir+"/part-00000", friction)
    p = lmgc2Interface.LmgcInterface()
else :
58
    p = scontact2.ParticleProblem(outputdir+"/part-00000")
dubois's avatar
maj  
dubois committed
59

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
60
deltamass = p.volume().sum()*(rhop-rho)
61
deltapinit = deltamass*-g/lxpart
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
62

dubois's avatar
maj  
dubois committed
63
64
print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
65
outf = 20 #000
dubois's avatar
maj  
dubois committed
66
67
68


# impermeability
69
70
strong_boundaries = [("Box",0,0),("Box",1,0),("Top",0,0.),("Top",1,0.),("Top",2,0.)]

Matthieu Constant's avatar
Matthieu Constant committed
71
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
72

dubois's avatar
maj  
dubois committed
73
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
74
fluid.export(outputdir,0,0)
dubois's avatar
maj  
dubois committed
75
76
77

tic = time.clock()
while t < tEnd :
78
    #if ii%100==0 and ii != 0:
Matthieu Constant's avatar
Matthieu Constant committed
79
     #  fluid.adapt_mesh(0.1,10,8e-4)
80
81
82
83
84
85
86
    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(10*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)  
dubois's avatar
maj  
dubois committed
87
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
88
    fluid.implicit_euler(dt)
dubois's avatar
maj  
dubois committed
89
90
91
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
92
93
        p.write(outputdir+"/part-%05d" % ioutput)
        fluid.export(outputdir, t, ioutput)
dubois's avatar
maj  
dubois committed
94
    ii += 1
95
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
dubois's avatar
maj  
dubois committed
96