avalanch2fluidsnofriction.py 4.71 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
Matthieu Constant's avatar
Matthieu Constant committed
23
24
25
26
27
28
29

# TESTCASE DESCRIPTION
# 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.

30
from migflow import fluid as fluid
Matthieu Constant's avatar
Matthieu Constant committed
31
from migflow import scontact
32
33
34
35
36
37
import numpy as np
import os
import time
import shutil
import random

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

Matthieu Constant's avatar
Matthieu Constant committed
42
43
44
45
46
47
48
49
50
# Physical parameters
g = -9.81                                        # gravity
rhop = 1500                                      # particles 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
51

Matthieu Constant's avatar
Matthieu Constant committed
52
53
54
55
# Numerical parameters
outf = 10                                        # number of iterations between output files
dt = 1e-4                                        # time step
tEnd = 1                                         # final time
56

Matthieu Constant's avatar
Matthieu Constant committed
57
58
59
60
61
62
# 
# PARTICLE PROBLEM
# 
# Particles structure builder to load file created with dep.py
p1 = scontact.ParticleProblem(2)
# Use particles deposit computed by the dep.py file of depot folder
63
p1.read_vtk("../depot/output",25)
Matthieu Constant's avatar
Matthieu Constant committed
64
65
# New particles structure builder to set other boundary conditions than in dep.py
p = scontact.ParticleProblem(2)
66
67
68
69
70
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)

Matthieu Constant's avatar
Matthieu Constant committed
71
72
73
# Initial time and iteration
t = 0
ii = 0
74

Matthieu Constant's avatar
Matthieu Constant committed
75
76
77
78
79
80
# 
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[num*rhom,nuf*rhof],[rhom,rhof],coeff_stab=0.001)
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
Nathan Coppin's avatar
Nathan Coppin committed
81
82
83
84
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Left")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Right")
Matthieu Constant's avatar
Matthieu Constant committed
85
# Set location of the particles in the mesh and compute the porosity in each computation cell
86
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
87
88
89
90
91
# Access to fluid fields and node coordiantes
s = fluid.solution()
x = fluid.coordinates()
# Set initial condition for the concentration
lx = 0.1
92
93
94
95
96
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
97
tic = time.time()
98
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
Matthieu Constant's avatar
Matthieu Constant committed
99
100
101
102

#
# COMPUTATION LOOP
# 
103
while t < tEnd : 
Matthieu Constant's avatar
Matthieu Constant committed
104
    # Fluid solver
105
    fluid.implicit_euler(dt)
106

Matthieu Constant's avatar
Matthieu Constant committed
107
    # Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
108
109
    if (ii%15==0 and ii != 0):
       fluid.adapt_mesh(1e-2,3e-3,10000)
110

Matthieu Constant's avatar
Matthieu Constant committed
111
    # Computation of the new velocities
112
    forces = fluid.compute_node_force(dt)
Matthieu Constant's avatar
Matthieu Constant committed
113
114
115
116
    vn = p.velocity() + forces * dt / p.mass()
    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
    # Number of sub time step
    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
117
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
Matthieu Constant's avatar
Matthieu Constant committed
118
    # NLGS iterations
119
120
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  
121

Matthieu Constant's avatar
Matthieu Constant committed
122
    t += dt
123
    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
124

Matthieu Constant's avatar
Matthieu Constant committed
125
    # Output files writting
126
127
128
129
130
    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
131
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))