boycott.py 10 KB
Newer Older
Margaux Boxho's avatar
Margaux Boxho 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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 14 14:49:05 2019

@author: margauxboxho
"""

# 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
# Bidimensional particles sedimentation in fluid


from migflow import fluid
from migflow import scontact

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

45

Margaux Boxho's avatar
Margaux Boxho committed
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
def genInitialPosition(filename, r, H, ly, lx, rhop, compacity) :
    """Set all the particles centre positions and create the particles objects to add in the computing structure
    
    Keyword arguments:    
    filename -- name of the output file
    r -- max radius of the particles
    H -- domain height
    ly - particles area height
    lx -- particles area width
    rhop -- particles density        
    """
    # 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
63
64
65
66
67
68
    #x = np.arange(-lx/2+r, lx/2-r, 2*r)
    #y = np.arange(0-r, -ly+r, -2*r)
    #x, y = np.meshgrid(x, y)
    #x = x.flat
    #y = y.flat
    
Margaux Boxho's avatar
Margaux Boxho committed
69
    N = 4*compacity*(lx*ly)/(np.pi*r**2)
70
71
72
73
74
    e = 2*ly/(N)**(1./2.)

    # Definition of the points where the particles are located
    shift = 0.1; alpha = 30*np.pi/180; 
    lx = lx*np.cos(alpha)
Margaux Boxho's avatar
Margaux Boxho committed
75
    
Margaux Boxho's avatar
Margaux Boxho committed
76
    x = np.arange(lx/2+25*r, -lx/2-25*r, -e)
77
78
79
80
81
    y = np.arange(-3*r-shift, -ly+3*r-shift, -e)
    
    #shift = 0.26;
    #x = np.arange(lx/2+50*r, -lx/2-50*r, -e)
    #y = np.arange(2*H-3*r-shift, 2*H-ly+3*r-shift, -e)
Margaux Boxho's avatar
Margaux Boxho committed
82
    x, y = np.meshgrid(x, y)
83
84
85
86
87
    
    
    x = x-y*np.tan(alpha); 
    y = y; 
    
Margaux Boxho's avatar
Margaux Boxho committed
88
89
90
    x = x.flat
    y = y.flat
    
91
    
Margaux Boxho's avatar
Margaux Boxho committed
92
    # Add a grain at each centre position
93
94
95
    for i in range(len(x)) :
        #p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
        p.add_particle((x[i]+2*random.uniform(-e/2+r,e/2-r), y[i]+2*random.uniform(-e/2+r,e/2-r)), r, r**2 * np.pi * rhop)
Margaux Boxho's avatar
Margaux Boxho committed
96
    p.write_vtk(filename,0,0)
97
    
Margaux Boxho's avatar
Margaux Boxho committed
98

99
100
filename = "/Volumes/Margaux/Memoire_Code/testcases/Boycott/"
outputdir = filename + "output"
Margaux Boxho's avatar
Margaux Boxho committed
101
102
103
104
105
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

# Physical parameters
g = -9.81                                       # gravity
106
107
r = 1e-3                                        # particles radius os sand seed
rhop = 1600                                     # particles density as sand
Margaux Boxho's avatar
Margaux Boxho committed
108
109
rho = 1000                                      # fluid density
nu = 1e-6                                       # kinematic viscosity
110
mu = rho*nu;        
111
compacity = 0.2
Margaux Boxho's avatar
Margaux Boxho committed
112
113

# Numerical parameters
Margaux Boxho's avatar
Margaux Boxho committed
114
outf = 5                                       # number of iterations between output files
115
116
dt = 0.5e-3                                    # time step
tEnd = 3                                       # final time
Margaux Boxho's avatar
Margaux Boxho committed
117
118

# Geometrical parameters
Margaux Boxho's avatar
Margaux Boxho committed
119
ly = 0.8                                        # particles area height
120
lx = 0.1                                        # particles area widht
121
H = 0.6                                         # domain height
Margaux Boxho's avatar
Margaux Boxho committed
122
123
124
125
126

#
# PARTICLE PROBLEM
#
# Initialise particles
127
#genInitialPosition(outputdir, r, H, ly, lx, rhop, compacity)
Margaux Boxho's avatar
Margaux Boxho committed
128
p = scontact.ParticleProblem(2)
129
p.read_vtk("/Volumes/Margaux/Memoire_Code/testcases/Boycott/output_no_slip_wall_Ref",0)
Margaux Boxho's avatar
Margaux Boxho committed
130
131
132
133

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

134
135
136
137
U_0 = 2*9.81*(rhop-rho)*r**2/(9*mu); 
phi = compacity
f_phi = ((1-phi)**2)/((1+phi**(1/3))*np.exp(5/3*phi/(1-phi))); 

138
139
U_0 = 0.8633; 
P_0 = 255;
140
141
142
Re = rho*U_0*r/mu; 
N_0 = len(p.position());
e_target = 1e-7
143
print(" ======= The initial number of grains is %.5f [-]=======" %N_0)
144
145
146
147
print(" ======= The initial velocity is %.5f [m/s]=======" %U_0)
print(" ======= The initial pressure is %.5f [Pa]=======" %P_0)
print(" ======= The Reynolds number is %.3f [Pa]=======" %Re)

Margaux Boxho's avatar
Margaux Boxho committed
148
149
150
151
# Initial time and iteration
t         =  0
ii        =  0
jj        =  0
152
153
delta     = 10

Margaux Boxho's avatar
Margaux Boxho committed
154
155
156
#
# FLUID PROBLEM
#
157
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho], petsc_solver_type="-pc_type lu -ksp_max_it 50")
Margaux Boxho's avatar
Margaux Boxho committed
158
159
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
160
161
162
fluid.set_wall_boundary("Bottom", velocity=[0,0])
fluid.set_wall_boundary("Lateral", velocity=[0,0])
fluid.set_wall_boundary("Top",pressure=0, velocity=[0,0])
Margaux Boxho's avatar
Margaux Boxho committed
163
# Set location of the particles in the mesh and compute the porosity in each computation cell
164
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
Margaux Boxho's avatar
Margaux Boxho committed
165
166
167
fluid.export_vtk(outputdir,0,0)

tic = time.time()
168
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
Margaux Boxho's avatar
Margaux Boxho committed
169
170
171
172
  
#
# COMPUTATION LOOP
#
173
while t < tEnd :
Margaux Boxho's avatar
Margaux Boxho committed
174
175
176
177
178
    if (ii==jj*delta):
        print("Save state at %i" %(ii))
        fluid.export_vtk(outputdir,0,0)
        p.write_vtk(outputdir,0,0)
        
179
180
181
    # Fluid solver
    fluid.implicit_euler(dt)
    
Margaux Boxho's avatar
Margaux Boxho committed
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
    if (ii==0):
        fluid.adapt_mesh(60*r, 5*r, 6602*2, e_target, 1, U_0, P_0)
        
    if (ii%(delta*(jj+1))!=0): 
        fluid.adapt_mesh(60*r, 5*r, 6602*2, e_target, 0, U_0, P_0) 
    
    # Adaptation of the mesh and loading the past save state
    if (ii%(delta*(jj+1))==0 and ii != 0):  
        print("Adapt mesh at %i" %(ii))
        # Adapt the mesh by creating the lc.pos file
        fluid.adapt_gen_mesh(60*r, 5*r, 6602*2, e_target, U_0, P_0)
        ii = ii-delta
        t = t-dt*delta
        jj += 1
        print("ii=%i, t=%i, jj=%i" %(ii,t,jj))
        # load the save past state 
198
199
        folder = outputdir+"/fluid_00000.vtu"
        fluid.import_vtk(folder)
Margaux Boxho's avatar
Margaux Boxho committed
200
201
202
203
204
205
206
207
        p.read_vtk(outputdir,0)
        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
        # create the .msh on the adapt/lc.pos
        fluid.adapt_load_mesh("adapt/mesh.msh")
        #fluid.implicit_euler(dt)
        print("Restart computation %i iterations back at %i" %(delta, ii)) 
        # Initialize the min_h_target attribut of the Mesh structure
        fluid.adapt_mesh(60*r, 5*r, 6602*2, e_target, 1, U_0, P_0) 
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
        
    # Computation of the new velocities
    forces = fluid.compute_node_force(dt)
    vn = p.velocity() + forces*dt/p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    # Number of contact sub time step
    nsub = max(1, int(np.ceil((vmax*dt*4)/min(p.r()))))
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
    # NLGS iterations
    for i in range(nsub) :
        tol = 1e-8                                #numerical tolerance for particles intersection
        p.iterate(dt/nsub, forces, tol)

    # Localisation of the particles in the fluid
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
    t += dt

    # Output files writting
    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))

"""
Margaux Boxho's avatar
Margaux Boxho committed
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296

while t < tEnd :
    # Save the state of the fluid at certain iteration 
    if (ii==jj*delta):
        print("Save state at %i" %(ii))
        fluid.export_vtk(outputdir,0,0)
        p.write_vtk(outputdir,0,0)
        
    # Fluid solver
    fluid.implicit_euler(dt)
    
    if (ii==0):
        fluid.adapt_mesh(500*r, 8*r, 6602*2, 0, 1)
        
    if (ii%(delta*(jj+1))!=0): 
        fluid.adapt_mesh(500*r, 8*r, 6602*2, 0, 0)
    
    # Adaptation of the mesh and loading the past save state
    if (ii%(delta*(jj+1))==0 and ii != 0):  
        print("Adapt mesh at %i" %(ii))
        #fluid.adapt_gen_mesh(500*r, 8*r, 6602*2, 1e-6, "adapt/lc.pos")
        # Adapt the mesh by creating the lc.pos file
        fluid.adapt_mesh(500*r, 8*r, 6602*2, 1, 0)
        ii = ii-delta
        t = t-dt*delta
        jj += 1
        print("ii=%i, t=%i, jj=%i" %(ii,t,jj))
        # load the save past state 
        fluid.import_vtk("output/fluid_00000.vtu")
        p.read_vtk(outputdir,0)
        fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
        # create the .msh on the adapt/lc.pos
        fluid.adapt_load_mesh("adapt/mesh.msh")
        #fluid.implicit_euler(dt)
        print("Restart computation %i iterations back at %i" %(delta, ii)) 
        # Initialize the min_h_target attribut of the Mesh structure
        fluid.adapt_mesh(500*r, 8*r, 6602*2, 0, 1)
    
    # Computation of the new velocities
    forces = fluid.compute_node_force(dt)
    # print(forces)
    vn = p.velocity() + forces*dt/p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    # Number of contact sub time step
    nsub = max(1, int(np.ceil((vmax*dt*4)/min(p.r()))))
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
    # NLGS iterations
    for i in range(nsub) :
        tol = 1e-8                                #numerical tolerance for particles intersection
        p.iterate(dt/nsub, forces, tol)

    # Localisation of the particles in the fluid
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
    t += dt

    # Output files writting
    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))

297
"""