depot-fluid.py 2.89 KB
Newer Older
dubois's avatar
maj  
dubois committed
1
#!/usr/bin/env python
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
2
use_lmgc=True
dubois's avatar
maj  
dubois committed
3
4
5
import pyfefp
from pyfefp.legendre import *
from pyfefp import  scontact2Interface
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
6
7
if use_lmgc :
    from pyfefp import  lmgc2Interface
dubois's avatar
maj  
dubois committed
8
9
10
11
12
13
14
import scontact2
import numpy as np
from pyfefp import mesh
import os
import time
import random

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


Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
27
outputdir = "output-initpressure"
dubois's avatar
maj  
dubois committed
28
29
30
31
32
33
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
34
r=1e-4*3
dubois's avatar
maj  
dubois committed
35
36
37
38
39
40
p = scontact2.ParticleProblem()

#physical parameters
g =  [0, -9.81]
rho = 1200
rhop = 1500 #2400
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
41
nu = 1e-6
dubois's avatar
maj  
dubois committed
42
43
44
tEnd = 20

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
47
c =  20 #h/dt* 7./50.
48
dt = 0.25e-4 #0.5*r/V
dubois's avatar
maj  
dubois committed
49
50
51
52
epsilon = 1e-5 #2e-2*c*h # ?? not sure ??

print("c = %g, eps = %g" % (c, epsilon))

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
53
y0part = 0.071
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
54
55
56
lypart = 0.025
lxpart = 0.015
genInitialPosition(outputdir + "/part-00000", r, 0., y0part, lxpart, lypart, rhop)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
57
58
59
60
61
62
if use_lmgc :
    friction=0.1
    lmgc2Interface.scontactToLmgc90(outputdir+"/part-00000", friction)
    p = lmgc2Interface.LmgcInterface()
else :
    p = scontact2Interface.ScontactInterface(outputdir+"/part-00000")
dubois's avatar
maj  
dubois committed
63

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
64
65
66
deltamass = p.volume().sum()*(rhop-rho)
deltapinit = deltamass*-g[1]/lxpart

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
71
fluid = pyfefp.FluidSolver("mesh.msh", dt, rhop=rhop, rho=rho, g = g, mu=nu*rho, c=c*c, epsilon = epsilon,timeScheme="Explicit", dragFactor=0.3, TypeEMu=TriangleP1, TypeEMp=TriangleP1)
dubois's avatar
maj  
dubois committed
72
73
74
75
76
77
78
79

# impermeability
fluid.add_boundary_condition("Top", 1, 0.)
fluid.add_boundary_condition("Top", 0, 0.)
fluid.add_boundary_condition("Box", 1, 0.)
fluid.add_boundary_condition("Box", 0, 0.)

fluid.complete()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
80
81
82
83
84
#y = fluid.mesh_position()[...,1].transpose()
#pinit = (y0part+lypart-y)*deltapinit/lypart
#pinit[pinit<0] = 0
#pinit[pinit>deltapinit] = deltapinit
#fluid._dof.setField(2, pinit, fluid.solution())
dubois's avatar
maj  
dubois committed
85
86
87
88
89
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

tic = time.clock()
while t < tEnd :
    forces = fluid.solve() + g * p.mass()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
90
    NS = 5 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
91
92
    if ii%NS == 0 :
        p.iterate(dt*NS, forces)
dubois's avatar
maj  
dubois committed
93
94
95
96
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
97
        #p.write_vtk(outputdir, ioutput, t)
dubois's avatar
maj  
dubois committed
98
        p.write(outputdir, ioutput, t)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
99
        #p.write_boundary_vtk(outputdir, ioutput, t)
100
        fluid.write_solution(outputdir, ioutput, t, "vtk")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
101
        #fluid.write_solution(outputdir, ioutput, t, "vtk")
dubois's avatar
maj  
dubois committed
102
103
104
    ii += 1
    print("%i : %.2g/%.2g (cpu %.2g)" % (ii, t, tEnd, time.clock() - tic))