depot.py 4.27 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/>.

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
37

import numpy as np
import os
import time
import shutil

def genInitialPosition(filename, rhop) :
Matthieu Constant's avatar
Matthieu Constant committed
38
39
40
41
42
43
44
45
    """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        
    """
    # Particles structure builder
46
    p = scontact.ParticleProblem(3)
47
    
Matthieu Constant's avatar
Matthieu Constant committed
48
    # Load the mesh.msh file specifying physical boundaries names
49
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "X", "Z"])
Matthieu Constant's avatar
Matthieu Constant committed
50
    # Reading the file containing the radii
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
51
    myRadius = np.genfromtxt('radius.csv', delimiter=',')
Nathan Coppin's avatar
Nathan Coppin committed
52
    r2 = np.amax(myRadius)*8e-6
Matthieu Constant's avatar
Matthieu Constant committed
53
54
    i = 0
    # Positioning the particles on a regular grid
Nathan Coppin's avatar
Nathan Coppin committed
55
    for y in np.arange(0, -.05, -2*r2):
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
56
57
        for z in np.arange(.01675-r2, -.01675, -2*r2):
            for x in np.arange(0.014-r2, -.014+r2, -2*r2):
Nathan Coppin's avatar
Nathan Coppin committed
58
                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
59
                i += 1
Matthieu Constant's avatar
Matthieu Constant committed
60
    # Translation of the particles to the top of the box
Nathan Coppin's avatar
Nathan Coppin committed
61
    p.position()[:, 1] += .07/2 
Matthieu Constant's avatar
Matthieu Constant committed
62
    print("Number of particles=%d" % len(p.position()[:, 1]))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
63
64
65
    p.write_vtk(filename, 0, 0)


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

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

Matthieu Constant's avatar
Matthieu Constant committed
77
# Numerical parameters
78
outf = 50                                        # number of iterations between output files
Matthieu Constant's avatar
Matthieu Constant committed
79
dt = 1e-4                                        # time step
Matthieu Constant's avatar
Matthieu Constant committed
80
tEnd = 100                                       # final time
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
81

Matthieu Constant's avatar
Matthieu Constant committed
82
83
84
85
#
# PARTICLE PROBLEM
#
# Particles structure creation and setting particles properties
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
86
genInitialPosition(outputdir, rhop)
Nathan Coppin's avatar
Nathan Coppin committed
87

88
p = scontact.ParticleProblem(3)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
89
90
p.read_vtk(outputdir,0)

Matthieu Constant's avatar
Matthieu Constant committed
91
92
93
94
95
96
97
# Initial time and iteration
t         =  0
ii        =  0

#
# FLUID PROBLEM
#
98

99
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho])
Matthieu Constant's avatar
Matthieu Constant committed
100
#Set the mesh geometry for the fluid computation
101
102
103
104
105
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
106
#Set location of the particles in the mesh and compute the porosity in each computation cell
107
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces(),init=True)
Matthieu Constant's avatar
Matthieu Constant committed
108
109

tic = time.time()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
110
fluid.export_vtk(outputdir,0,0)
Matthieu Constant's avatar
Matthieu Constant committed
111

Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
112
G = np.zeros((p.n_particles(),p.dim())
113
G[:,1] = p.mass()[:,0]*g[1]
Matthieu Constant's avatar
Matthieu Constant committed
114
115
116
#
# COMPUTATION LOOP
#
Matthieu Constant's avatar
Matthieu Constant committed
117
while t < tEnd : 
118
    time_integration.iterate(fluid, p, dt,min_nsub=20, contact_tol = 1e-8,external_particles_forces=G)
Matthieu Constant's avatar
Matthieu Constant committed
119
120
    t += dt
    # Output files writing
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
121
122
123
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)
Nathan Coppin's avatar
Nathan Coppin committed
124
        #fluid.export_vtk(outputdir, t, ioutput)
Matthieu Constant's avatar
Matthieu Constant committed
125
    ii += 1
Matthieu Constant's avatar
Matthieu Constant committed
126
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))