depot.py 4.15 KB
Newer Older
Michel Henry's avatar
Michel Henry committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
#
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
# 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
# GNU Lesser General Public License for more details.
#
# 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,
# see <http://www.gnu.org/licenses/>.

#!/usr/bin/env python


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

from migflow import fluid
from migflow import scontact
from migflow import time_integration

import numpy as np
import os
import time
import shutil

def genInitialPosition(p, filename, rhop) :
    """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

    Keyword arguments:
    filename -- name of the output file
    rhop -- particles density
    """
    # Reading the file containing the radii
    myRadius = np.genfromtxt('radius.csv', delimiter=',')*2
    r2 = np.amax(myRadius)*8e-6
    i = 0
    # Positioning the particles on a regular grid
    for y in np.arange(0, -.05, -2*r2):
        for z in np.arange(.01675-r2, -.01675, -2*r2):
            for x in np.arange(0.014-r2, -.014+r2, -2*r2):
                p.add_particle((x, y, z), myRadius[i%500000]*8e-6, 4./3.* (myRadius[i%500000]*8e-6)**3 * np.pi * rhop);
                i += 1
    # Translation of the particles to the top of the box
    state = p.state()
    state.x[:, 1] += .07/2
    p.set_state(state)
    print("Number of particles=%d" % len(p.state().x[:, 1]))
    # p.write_vtk(filename, 0, 0)


# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

# Physical parameters
g = np.array([0,-9.81,0])                        # gravity
rho = 1000                                       # fluid density
rhop = 2640                                      # particles density
nu = 1e-6                                        # kinematic viscosity

# Numerical parameters
outf = 50                                        # number of iterations between output files
dt = 1e-4                                        # time step
tEnd = 100                                       # final time

#
# PARTICLE PROBLEM
#
# Particles structure creation and setting particles properties

p = scontact.ParticleProblem(3)
p.load_msh_boundaries("mesh.msh", ["Bottom", "Top", "X", "Z"])
genInitialPosition(p,outputdir, rhop)
print(p.n_particles())
p.write_vtk(outputdir,0,0)

# Initial time and iteration
t         =  0
ii        =  0

#
# FLUID PROBLEM
#

fluid = fluid.FluidProblem(3,g,[nu*rho],[rho])
#Set the mesh geometry for the fluid computation
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")
#Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces(),init=True)

tic = time.time()
fluid.export_vtk(outputdir,0,0)
G = np.zeros((p.n_particles(),p.dim()))
G[:,1] = p.mass()[:,0]*g[1]
#
# COMPUTATION LOOP
#
while t < tEnd :
    time_integration.iterate(fluid, p, dt,min_nsub=20, contact_tol = 1e-8,external_particles_forces=G)
    t += dt
    # Output files writing
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
        #fluid.export_vtk(outputdir, t, ioutput)
    ii += 1
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))