avalanch2fluids.py

# 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
# 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. 
# particles are initially placed using a random loose packing method of the lmgc90 software.
# The lmgc90 software is also used to solve contacts in this test cases.

from migflow import fluid
from migflow import scontact
from migflow import lmgc90Interface
from pylmgc90 import pre
import numpy as np
import os
import time
import shutil
import random

def genInitialPosition(filename, N, r, lx, ly, rhop) :
    """Create a container to initialize particles and write particles in an initial output file.

    Keyword arguments:
    filename -- directory to write the output file
    N -- maximum number of grains
    r -- mean particle radius
    lx -- particle area width
    ly -- particle aera height
    rhop -- particle density
    """
    # Particles structure builder
    p = scontact.ParticleProblem(2)
    # Load mesh.msh file specifying physical boundaries names
    p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Right", "Left"])
    # Radii defining the interval of the granulo
    Ra = 0.9*r
    Rb = 1.1*r
    # Random list of radii in the interval
    radii = pre.granulo_Random(N, Ra, Rb)
    # Nb of particles and positions in the geometry defined by arguments: [0,lx*4]*[0,ly/10]
    [nb_laid_particles, coors] = pre.depositInBox2D(radii,lx*4,ly/10)
    # Add a grain at each centre position
    for i in range(nb_laid_particles):
       p.add_particle(coors[2*i:2*i+2], radii[i], radii[i]**2 * np.pi * rhop)
    # Particles horizontal translation
    p.position()[:, 0] += 0.1

    # Second particles area
    radii = pre.granulo_Random(N, Ra, Rb)
    [nb_laid_particles, coors] = pre.depositInBox2D(radii,lx,ly)
    for i in range(nb_laid_particles):
       p.add_particle(coors[2*i:2*i+2], radii[i], radii[i]**2 * np.pi * rhop)
    p.write_vtk(filename,0,0)

# Define output directory
outputdir = "output2fluids"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

# Physical parameters
g = -9.81                                     # gravity
# Particles
r = 0.0005                                    # radius
lx = 0.1                                      # particles area width
ly = 0.18                                     # particles area height
N = 15000                                     # number of particles
rhop = 1100                                   # solid density
# Soft water
rhof = 1000                                   # soft water density
nuf = 1e-6                                    # soft water viscosity
# Mixture properties
rhom = 1050                                   # sea water density
num = 1e-6                                    # sea water kinematic viscosity

# Numerical parameters
outf = 10                                     #number of iterations between output files
dt = 1e-4                                     #time step
tEnd = 10                                     #final time

# 
# PARTICLE PROBLEM
# 
# Initialise particles
genInitialPosition(outputdir, N, r, lx, ly, rhop)
friction = 0.3                                #friction coefficient
lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
p = lmgc90Interface.ParticleProblem(2)

# Initial time and iteration
t = 0
ii = 0

# 
# FLUID PROBLEM
# 
fluid = fluid.FluidProblem(2,g,[num*rhom,nuf*rhof],[rhom,rhof])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Left")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Right")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#Access to fluid fields and node coordiantes
s         =  fluid.solution()
x         =  fluid.coordinates()
#Set initial condition for the concentration
for i in range(len(x[:,0])):
    if (x[i,0])<lx:
        s[i,3] = 1
fluid.export_vtk(outputdir,0,0)

tic       =  time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())

# 
# COMPUTATION LOOP
# 
while t < tEnd : 
    # Fluid solver
    fluid.implicit_euler(dt)

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

    # 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()))
    # Contact solver
    for i in range(nsub) :
        p.iterate(dt/nsub, forces)  

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

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