depot.py 4.06 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
Michel Henry's avatar
Michel Henry 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
# 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
Jonathan Lambrechts's avatar
Jonathan Lambrechts 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
#
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
13
14
15
# 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.
Michel Henry's avatar
Michel Henry committed
17
#
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,
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
21
# see <http://www.gnu.org/licenses/>.

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
22
#!/usr/bin/env python
Matthieu Constant's avatar
Matthieu Constant committed
23
24
25
26
27
28
29


# TESTCASE DESCRIPTION
# Deposit of 3D particles in a box full of fluid.

from migflow import fluid
from migflow import scontact
30
from migflow import time_integration
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
31
32
33
34
35
36

import numpy as np
import os
import time
import shutil

Michel Henry's avatar
Michel Henry committed
37
def genInitialPosition(p, filename, rhop) :
Matthieu Constant's avatar
Matthieu Constant committed
38
39
    """Set all the particles centre positions and create the particles objects to add in the computing structure.
    For this testcase, radius are contained in a file radius.csv created with a special distribution
Michel Henry's avatar
Michel Henry committed
40

Matthieu Constant's avatar
Matthieu Constant committed
41
42
    Keyword arguments:
    filename -- name of the output file
Michel Henry's avatar
Michel Henry committed
43
    rhop -- particles density
Matthieu Constant's avatar
Matthieu Constant committed
44
45
    """
    # Reading the file containing the radii
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
46
    myRadius = np.genfromtxt('radius.csv', delimiter=',')
Nathan Coppin's avatar
Nathan Coppin committed
47
    r2 = np.amax(myRadius)*8e-6
Matthieu Constant's avatar
Matthieu Constant committed
48
49
    i = 0
    # Positioning the particles on a regular grid
Nathan Coppin's avatar
Nathan Coppin committed
50
    for y in np.arange(0, -.05, -2*r2):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
51
52
        for z in np.arange(.01675-r2, -.01675, -2*r2):
            for x in np.arange(0.014-r2, -.014+r2, -2*r2):
Michel Henry's avatar
Michel Henry committed
53
                p.add_particle((x, y, z), myRadius[i%500000]*8e-6, 4./3.* (myRadius[i%500000]*8e-6)**3 * np.pi * rhop);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
54
                i += 1
Matthieu Constant's avatar
Matthieu Constant committed
55
    # Translation of the particles to the top of the box
56
57
    p.position()[:,1] += .07/2
    print("Number of particles=%d" % p.position().shape[0])
Michel Henry's avatar
Michel Henry committed
58
    # p.write_vtk(filename, 0, 0)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
59
60


Matthieu Constant's avatar
Matthieu Constant committed
61
# Define output directory
62
outputdir = "output"
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
63
64
65
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

Matthieu Constant's avatar
Matthieu Constant committed
66
# Physical parameters
67
g = np.array([0,-9.81,0])                        # gravity
Matthieu Constant's avatar
Matthieu Constant committed
68
69
70
rho = 1000                                       # fluid density
rhop = 2640                                      # particles density
nu = 1e-6                                        # kinematic viscosity
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
71

Matthieu Constant's avatar
Matthieu Constant committed
72
# Numerical parameters
73
outf = 50                                        # number of iterations between output files
Matthieu Constant's avatar
Matthieu Constant committed
74
dt = 1e-4                                        # time step
Matthieu Constant's avatar
Matthieu Constant committed
75
tEnd = 100                                       # final time
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
76

Matthieu Constant's avatar
Matthieu Constant committed
77
78
79
80
#
# PARTICLE PROBLEM
#
# Particles structure creation and setting particles properties
Nathan Coppin's avatar
Nathan Coppin committed
81

82
p = scontact.ParticleProblem(3)
Michel Henry's avatar
Michel Henry committed
83
84
85
86
p.load_msh_boundaries("mesh.msh", ["Bottom", "Top", "X", "Z"])
genInitialPosition(p,outputdir, rhop)
print(p.n_particles())
p.write_vtk(outputdir,0,0)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
87

Matthieu Constant's avatar
Matthieu Constant committed
88
89
90
91
92
93
94
# Initial time and iteration
t         =  0
ii        =  0

#
# FLUID PROBLEM
#
95

96
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho])
Matthieu Constant's avatar
Matthieu Constant committed
97
#Set the mesh geometry for the fluid computation
98
99
100
101
102
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("X")
fluid.set_wall_boundary("Z")
Matthieu Constant's avatar
Matthieu Constant committed
103
#Set location of the particles in the mesh and compute the porosity in each computation cell
104
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
Matthieu Constant's avatar
Matthieu Constant committed
105
106

tic = time.time()
Nathan Coppin's avatar
Nathan Coppin committed
107
fluid.write_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
108
109
110
#
# COMPUTATION LOOP
#
Michel Henry's avatar
Michel Henry committed
111
while t < tEnd :
Nathan Coppin's avatar
Nathan Coppin committed
112
    time_integration.iterate(fluid, p, dt,min_nsub=20, contact_tol = 1e-8,external_particles_forces=g*p.mass())
Matthieu Constant's avatar
Matthieu Constant committed
113
114
    t += dt
    # Output files writing
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
115
116
117
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
Nathan Coppin's avatar
Nathan Coppin committed
118
        fluid.write_vtk(outputdir, ioutput, t)
Matthieu Constant's avatar
Matthieu Constant committed
119
    ii += 1
Matthieu Constant's avatar
Matthieu Constant committed
120
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))