depot.py 4.36 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
Matthieu Constant's avatar
Matthieu Constant committed
2
3
# <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

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

Michel Henry's avatar
Michel Henry committed
38
def genInitialPosition(filename, r, H, ly, lx, rhop) :
Matthieu Constant's avatar
Matthieu Constant committed
39
    """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
    """
    # Particles structure builder
Michel Henry's avatar
Michel Henry committed
50
    p = scontact.ParticleProblem(2)
Matthieu Constant's avatar
Matthieu Constant committed
51
52
53
54
    # 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
Michel Henry's avatar
Michel Henry committed
66
outputdir = "output_incompressible"
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
Nathan Coppin's avatar
Nathan Coppin committed
72
r = 1.5e-3                                        # particles radius
Matthieu Constant's avatar
Matthieu Constant committed
73
74
75
76
77
rhop = 1500                                     # particles density
rho = 1000                                      # fluid density
nu = 1e-6                                       # kinematic viscosity

# Numerical parameters
Michel Henry's avatar
Michel Henry committed
78
outf = 2                                        # number of iterations between output files
Nathan Coppin's avatar
Nathan Coppin committed
79
dt = 1e-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

#
# PARTICLE PROBLEM
#
# Initialise particles
Michel Henry's avatar
Michel Henry committed
91
genInitialPosition(outputdir, r, H, ly, lx, rhop)
Matthieu Constant's avatar
Matthieu Constant committed
92
p = scontact.ParticleProblem(2)
Michel Henry's avatar
Michel Henry committed
93
p.read_vtk(outputdir,0)
Matthieu Constant's avatar
Matthieu Constant committed
94
95
96
97
98
99
100
101
102
103
104

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
#
Michel Henry's avatar
Michel Henry committed
105
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=1,solver="petsc4py",solver_options="-pc_type lu",usolid=True)
Matthieu Constant's avatar
Matthieu Constant committed
106
107
# 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.position(), p.velocity(), p.contact_forces())
Nathan Coppin's avatar
Nathan Coppin committed
113
fluid.write_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
118
119
120

tic = time.time()
#
# COMPUTATION LOOP
#
while t < tEnd :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
121
    time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.mass(), use_predictor_corrector=True)
Matthieu Constant's avatar
Matthieu Constant committed
122
123
124
125
126
    t += dt

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