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);
Matthieu Constant's avatar
Matthieu Constant committed
21
    p.write_vtk(filename,0,0)
dubois's avatar
maj  
dubois committed
22
23


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
lxpart = 0.015
Matthieu Constant's avatar
Matthieu Constant committed
52
genInitialPosition(outputdir, 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 :
Matthieu Constant's avatar
Matthieu Constant committed
58
59
    p = scontact2.ParticleProblem()
    p.read_vtk(outputdir,0)
dubois's avatar
maj  
dubois committed
60

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

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


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

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

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

tic = time.clock()
Matthieu Constant's avatar
Matthieu Constant committed
78
79
80
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd : 
    fluid.implicit_euler(dt)
81
82
83
    forces = fluid.compute_node_force(dt)
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
84
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
85
86
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
87
        p.iterate(dt/nsub, forces,1e-6)  
dubois's avatar
maj  
dubois committed
88
89
90
91
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Matthieu Constant's avatar
Matthieu Constant committed
92
93
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(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