sablier.py 4.53 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2018>
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
4
# 	
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
5
# List of the contributors to the development of MigFlow: see AUTHORS file.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
6
7
# Description and complete License: see LICENSE file.
# 	
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
8
# This program (MigFlow) is free software: 
9
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
10
11
12
13
14
15
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16
# GNU Lesser General Public License for more details.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
17
# 
18
19
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files).  If not, 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
21
# see <http://www.gnu.org/licenses/>.

dubois's avatar
maj  
dubois committed
22
#!/usr/bin/env python
Matthieu Constant's avatar
Matthieu Constant committed
23
24
25
26
27

# TESTCASE DESCRIPTION
# Falling particles in a sandglass
from migflow import fluid
from migflow import scontact
Frédéric Dubois's avatar
Frédéric Dubois committed
28
from migflow import lmgc90Interface
29

dubois's avatar
maj  
dubois committed
30
31
32
import numpy as np
import os
import time
33
import shutil
dubois's avatar
maj  
dubois committed
34
import random
Matthieu Constant's avatar
Matthieu Constant committed
35
36

# grains creation and placement + physical boundaries definition
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
37
def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
Matthieu Constant's avatar
Matthieu Constant committed
38
    p = scontact.ParticleProblem(2)
Matthieu Constant's avatar
Matthieu Constant committed
39
    #Loading of the mesh.msh file specifying physical boundaries name
40
    p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
Matthieu Constant's avatar
Matthieu Constant committed
41
    #Definition of the points where the grains are located
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
42
43
    x = np.arange(ox+r, ox+lx, 2*r)
    y = np.arange(oy+r, oy+ly, 2*r)
dubois's avatar
maj  
dubois committed
44
45
    for i in range(x.shape[0]) :
       for j in range(y.shape[0]) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
46
          rr=r*(0.95+random.random()*0.05)
Matthieu Constant's avatar
Matthieu Constant committed
47
          #Addition of an particle object at each point
48
          p.add_particle((x[i], y[j]), rr, rr**2 * np.pi * rhop);
Matthieu Constant's avatar
Matthieu Constant committed
49
    p.write_vtk(filename,0,0)
dubois's avatar
maj  
dubois committed
50
51


52
outputdir = "output"
dubois's avatar
maj  
dubois committed
53
54
55
56
57
58
59
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0

#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
60
61
62
63
64
65
66
67
68
g =  -9.81                                      # gravity
rho = 1200                                      # fluid density
rhop = 1800                                     # grains density
nu = 1e-6                                       # kinematic viscosity
tEnd = 7                                        # final time
r=1e-4                                          # grains radius
y0part = 40e-3                                  # vertical center of the grains
lypart = 20e-3                                  # semi-height of the grains
lxpart = 0.015                                  # semi-width of the grains
dubois's avatar
maj  
dubois committed
69
70

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
71
dt = 1e-3                                       # time step
dubois's avatar
maj  
dubois committed
72

Matthieu Constant's avatar
Matthieu Constant committed
73
#Enable the use of lmgc90
Frédéric Dubois's avatar
Frédéric Dubois committed
74
use_lmgc90 = True
Matthieu Constant's avatar
Matthieu Constant committed
75
76

#Object particles creation
Matthieu Constant's avatar
Matthieu Constant committed
77
genInitialPosition(outputdir, r, 0., y0part, lxpart, lypart, rhop)
Frédéric Dubois's avatar
Frédéric Dubois committed
78
if use_lmgc90 :
Matthieu Constant's avatar
Matthieu Constant committed
79
    friction=0.1                                # friction coefficient
Frédéric Dubois's avatar
Frédéric Dubois committed
80
81
    lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
    p = lmgc90Interface.ParticleProblem(2)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
82
else :
Matthieu Constant's avatar
Matthieu Constant committed
83
    p = scontact.ParticleProblem(2)
Matthieu Constant's avatar
Matthieu Constant committed
84
    p.read_vtk(outputdir,0)
dubois's avatar
maj  
dubois committed
85
86

print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
Matthieu Constant's avatar
Matthieu Constant committed
87
outf = 20                                       # number of iterations between output files
dubois's avatar
maj  
dubois committed
88
89


Matthieu Constant's avatar
Matthieu Constant committed
90
91
92
93
94
95
96
fluid = fluid.FluidProblem(2,g,nu*rho,rho)
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Box",0,0)
fluid.set_strong_boundary("Box",1,0)
fluid.set_strong_boundary("Box",0,0)
fluid.set_strong_boundary("Box",1,0)
fluid.set_strong_boundary("Box",2,0)
dubois's avatar
maj  
dubois committed
97
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
98
fluid.export_vtk(outputdir,0,0)
dubois's avatar
maj  
dubois committed
99

Matthieu Constant's avatar
Matthieu Constant committed
100
tic = time.time()
Matthieu Constant's avatar
Matthieu Constant committed
101
#Computation loop
Matthieu Constant's avatar
Matthieu Constant committed
102
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
103
104
while t < tEnd :
    #Fluid solver
Matthieu Constant's avatar
Matthieu Constant committed
105
    fluid.implicit_euler(dt)
106
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
107
    #Computation of the new velocities
108
109
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
Matthieu Constant's avatar
Matthieu Constant committed
110
    #number of sub time step
Matthieu Constant's avatar
Matthieu Constant committed
111
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
112
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
113
    #Contact solver
114
    for i in range(nsub) :
Matthieu Constant's avatar
Matthieu Constant committed
115
        tol = 1e-8
Matthieu Constant's avatar
Matthieu Constant committed
116
117
        p.iterate(dt/nsub, forces, tol)
    #Localisation of the grains in the fluid
dubois's avatar
maj  
dubois committed
118
119
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt
Matthieu Constant's avatar
Matthieu Constant committed
120
    #Output files writting
dubois's avatar
maj  
dubois committed
121
122
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Matthieu Constant's avatar
Matthieu Constant committed
123
124
        p.write_vtk(outputdir, ioutput, t)
        fluid.export_vtk(outputdir, t, ioutput)
dubois's avatar
maj  
dubois committed
125
    ii += 1
Frédéric Dubois's avatar
Frédéric Dubois committed
126
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))