avalanch2fluidsnofriction.py 5.95 KB
Newer Older
1
# MigFlow - Copyright (C) <2010-2018>
2
3
4
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
5
# List of the contributors to the development of MigFlow: see AUTHORS file.
6
7
# Description and complete License: see LICENSE file.
# 	
8
# This program (MigFlow) is free software: 
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# 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
23
24
25
26
27
'''
Collapsing of a grain column in fluid. The grain column is initially mixed in sea water while the rest of the domain contain soft water. 
This test case cannot be runned on its own, it requires the output files of the depot test case (../depot/dep.py) to create a compact column as initial situation
Positions, masses and radii of graisn are setted in the dep.py file.
'''
28
29
30
31
32
33
34
35
36
from migflow import fluid as fluid
from migflow import scontact2
from pylmgc90 import pre
import numpy as np
import os
import time
import shutil
import random

37
outputdir = "output"
38
39
40
41
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

#physical parameters
42
43
44
45
46
47
48
49
g         = -9.81                                      #gravity
rhop      =  1500                                      #grains density
##Soft water
rhof      =  1000                                      #sof water density
nuf       =  1e-6                                      #soft water kinematic viscosity
##sea water
rhom      =  1050                                      #sea water density
num       =  1e-6                                      #sea water kinematic viscosity
50
51

#numerical parameters
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
outf      =  10                                        #number of iterations between output files
dt        =  1e-3                                      #time step
tEnd      =  1                                         #final time

#Grains structure builder to load file created with dep.py
p1        = scontact2.ParticleProblem()
#Use grains deposit computed by the dep.py file of depot folder
p1.read_vtk("../depot/output",25)
#New Grains structure builder to set other boundary conditions than in dep.py
p         = scontact2.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Left", "Right"])
for i in range(len(p1.position()[:,0])):
    p.add_particle((p1.position()[i,0], p1.position()[i,1]), p1.r()[i], p1.r()[i]**2 * np.pi * rhop)
p.write_vtk(outputdir, 0, 0)
#Initial time and iteration
t         =  0
ii        =  0

'''
fluid_problem is the function builder of the fluid class
Args:     -g             : gravity
          -viscosity     : list of fluid phases viscosities
          -density       : list of fluid phases densities
          -coeffStab     : stabilisation parametre to limit overshoots in concentration field (optional)
'''
fluid     =  fluid.fluid_problem(g,[num*rhom,nuf*rhof],[rhom,rhof],coeffStab=0.001)
#Set the mesh geometry for the fluid computation
79
80
fluid.load_msh("mesh.msh")

81
82
83
84
85
86
'''
set_strong_boundary is the function setting the strong boundaries (=constraints) for the fluid problem
Args:     -Tag           : physical tag (set in the mesh.geo file) of the boundary on which you want to add a constraint
          -Field         : O: x-velocity; 1: y-velocity; 2: pressure
          -Value         : value assigned to the specified field on the specified boundary
'''
87
88
89
90
91
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Right",0,0)
fluid.set_strong_boundary("Left",0,0)
92
93
94
95
96
'''
set_weak_boundary is the function setting the weak boundaries (=normal fluxes) for the fluid problem
Args:     -Tag           : physical tag (set in the mesh.geo file) of the boundary on which you want to add a constraint
          -Bnd Type      : boundaries can be of type: wall, inflow, outflow
'''
97
98
99
100
101
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Left","Wall")
fluid.set_weak_boundary("Top","Wall")
fluid.set_weak_boundary("Right","Wall")

102
#Set location of the grains in the mesh and compute the porosity in each computation cell
103
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
104
105
106
107
108
#Access to fluid fields and node coordiantes
s         =  fluid.solution()
x         =  fluid.coordinates()
#Set initial condition for the concentration
lx        =  0.1
109
110
111
112
113
for i in range(len(x[:,0])):
    if (x[i,0])<lx:
        s[i,3] = 1
fluid.export_vtk(outputdir,0,0)

Matthieu Constant's avatar
Matthieu Constant committed
114
tic       =  time.time()
115
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
116
#Computation loop
117
118
119
while t < tEnd : 
    #Fluid solver
    fluid.implicit_euler(dt)
120
121

    #Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
122
123
    if (ii%15==0 and ii != 0):
       fluid.adapt_mesh(1e-2,3e-3,10000)
124

125
    #Computation of the new velocities
126
127
128
    forces = fluid.compute_node_force(dt)
    vn     = p.velocity() + forces * dt / p.mass()
    vmax   = np.max(np.hypot(vn[:, 0], vn[:, 1]))
129
    #number of sub time step
130
    nsub   = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
131
132
133
134
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
    #Contact solver
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  
135
136

    t     += dt
137
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
138

139
140
141
142
143
144
    #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
Matthieu Constant's avatar
Matthieu Constant committed
145
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))