depot.py 4.29 KB
Newer Older
Matthieu Constant's avatar
Matthieu Constant committed
1
2
3
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
Michel Henry's avatar
Michel Henry committed
4
#
Matthieu Constant's avatar
Matthieu Constant committed
5
6
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
Michel Henry's avatar
Michel Henry committed
7
8
9
#
# This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
Matthieu Constant's avatar
Matthieu Constant committed
10
11
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
Michel Henry's avatar
Michel Henry committed
12
#
Matthieu Constant's avatar
Matthieu Constant committed
13
14
15
16
# 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
# GNU Lesser General Public License for more details.
Michel Henry's avatar
Michel Henry committed
17
#
Matthieu Constant's avatar
Matthieu Constant committed
18
# You should have received a copy of the GNU Lesser General Public License
Michel Henry's avatar
Michel Henry committed
19
# along with this program (see COPYING and COPYING.LESSER files).  If not,
Matthieu Constant's avatar
Matthieu Constant committed
20
21
22
23
24
25
26
27
28
29
# see <http://www.gnu.org/licenses/>.

#!/usr/bin/env python

# TESTCASE DESCRIPTION
# Bidimensional particles sedimentation in fluid


from migflow import fluid
from migflow import scontact
30
from migflow import time_integration
Matthieu Constant's avatar
Matthieu Constant committed
31
32
33
34
35
36
37
38
39

import numpy as np
import os
import time
import shutil
import random

def genInitialPosition(filename, r, H, ly, lx, rhop) :
    """Set all the particles centre positions and create the particles objects to add in the computing structure
Michel Henry's avatar
Michel Henry committed
40
41

    Keyword arguments:
Matthieu Constant's avatar
Matthieu Constant committed
42
43
44
45
46
    filename -- name of the output file
    r -- max radius of the particles
    H -- domain height
    ly - particles area height
    lx -- particles area width
Michel Henry's avatar
Michel Henry committed
47
    rhop -- particles density
Matthieu Constant's avatar
Matthieu Constant committed
48
49
50
51
52
53
54
    """
    # Particles structure builder
    p = scontact.ParticleProblem(2)
    # Load mesh.msh file specifying physical boundaries names
    p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom"])

    #Definition of the points where the particles are located
55
    x = np.arange(-lx/2+r-1e-8, lx/2-r+1e-8, 2*r)
Matthieu Constant's avatar
Matthieu Constant committed
56
57
58
59
60
61
62
63
64
65
    y = np.arange(H/2-r, H/2-ly+r, -2*r)
    x, y = np.meshgrid(x, y)
    x = x.flat
    y = y.flat
    # Add a grain at each centre position
    for i in range(len(x)) :
        p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
    p.write_vtk(filename,0,0)

# Define output directory
66
outputdir = "output"
Matthieu Constant's avatar
Matthieu Constant committed
67
68
69
70
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

# Physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
71
g = np.array([0,-9.81])                         # gravity
Matthieu Constant's avatar
Matthieu Constant committed
72
73
74
75
76
77
78
r = 1e-3                                        # particles radius
rhop = 1500                                     # particles density
rho = 1000                                      # fluid density
nu = 1e-6                                       # kinematic viscosity

# Numerical parameters
outf = 5                                        # number of iterations between output files
Matthieu Constant's avatar
Matthieu Constant committed
79
dt = 2.5e-3                                     # time step
Matthieu Constant's avatar
Matthieu Constant committed
80
81
82
tEnd = 100                                      # final time

# Geometrical parameters
Matthieu Constant's avatar
Matthieu Constant committed
83
84
ly = 1e-1                                       # particles area height
lx = 4e-1                                       # particles area widht
85
H = 0.6                                         # domain height
Matthieu Constant's avatar
Matthieu Constant committed
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, H, ly, lx, rhop)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)

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

# Initial time and iteration
t         =  0
ii        =  0

#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
Matthieu Constant's avatar
Matthieu Constant committed
108
109
110
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
fluid.set_wall_boundary("Top",pressure=0)
Matthieu Constant's avatar
Matthieu Constant committed
111
# Set location of the particles in the mesh and compute the porosity in each computation cell
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
112
fluid.set_particles(p.mass(), p.volume(), p.state(), p.contact_forces(),init=True)
Matthieu Constant's avatar
Matthieu Constant committed
113
fluid.export_vtk(outputdir,0,0)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
114
p.write_vtk(outputdir, 0, 0)
Matthieu Constant's avatar
Matthieu Constant committed
115
116
117

tic = time.time()

Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
118
G = np.zeros((p.n_particles(),2))
119
G[:,1] = p.mass()[:,0]*g[1]
Matthieu Constant's avatar
Matthieu Constant committed
120
121
122
123
#
# COMPUTATION LOOP
#
while t < tEnd :
124
    time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=G)
Matthieu Constant's avatar
Matthieu Constant committed
125
126
127
128
129
    t += dt

    # Output files writting
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
130
        p.write_vtk(outputdir, ioutput, t)
Matthieu Constant's avatar
Matthieu Constant committed
131
132
133
        fluid.export_vtk(outputdir, t, ioutput)
    ii += 1
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))