Download :download:`this example <../../examples/avalanche.py>`.

.. code-block:: python
  
  # MigFlow - Copyright (C) <2010-2026>
  # <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/>.
  

Bidimensional elongated particles avalanch into a air-water box.
================================================================
This example described how to initialise a particle problem with elongated bodies with prescribed orientation.
It also describes how to create a 2-fluids problem into MigFlow, air and water are considered in this example.

.. code-block:: python
  
  import os
  import gmsh
  import time
  import shutil
  import random
  import numpy as np
  from migflow import fluid
  from migflow import scontact
  from migflow import time_integration
  

First, let's create the output directory.
Define output directory

.. code-block:: python
  
  outputdir = "output-avalanche"
  shutil.rmtree(outputdir,ignore_errors=True)
  if not os.path.isdir(outputdir) :
      os.makedirs(outputdir)

Domain parameters :

.. code-block:: python
  
  height = 0.3                                                # domain height
  width  = 1.0                                                # domain width

Body description
----------------
We will see how to create a initial packing with a chosen bodies orientation
The body is described by its mask over a grid.

.. code-block:: python
  
  model = np.array([[1,1,1,1]])
  n_p_b = np.sum(model)                                       # number of particles by body
  angle = np.pi/4

Lets's define the body and particles properties :

.. code-block:: python
  
  g        = np.array([0, -9.81])                             # gravity
  rhop     = 1840                                             # particle density
  friction = 0.1                                              # friction coefficient
  r        = 2e-3                                             # maximal particle radius
  d        = 2*r                                              # maximal particle diameter
  b_size   = d*n_p_b                                          # maximal body size

Initial disk particles depot
----------------------------

First, we consider the falling of bidimensional particles into an open box.
Each particle has a diameter of the body size.

.. code-block:: python
  
  init_disk = scontact.ParticleProblem(2)
  init_disk.set_fixed_contact_geometry(0)
  init_disk.set_friction_coefficient(friction)
  

The open box is defined using three segments.

.. code-block:: python
  
  w  = 0.075                                                  # open-box width
  h  = 2.5                                                    # open-box height
  init_disk.add_boundary_segment([0, 0], [0, h], "left")
  init_disk.add_boundary_segment([w, h], [w, 0], "right")
  init_disk.add_boundary_segment([w, 0], [0, 0], "bottom")

The particles are then generated on a compact rectangular grid of size lx X ly :

.. code-block:: python
  
  x0     = np.array([w/2, h/2])
  lx     = .5*w
  ly     = .75*h
  x      = np.arange(-lx/2, lx/2, b_size) + x0[0]
  y      = np.arange(-ly/2, ly/2, b_size) + x0[1]
  x,y    = np.meshgrid(x,y)
  bodies = np.asarray([], int)
  for xi,yi in zip(x.flat, y.flat):
      ri = .5*np.random.uniform(0.8,1)*b_size
      body = init_disk.add_particle((xi,yi), ri, ri**2*np.pi*rhop)
      bodies = np.append(bodies, body)

Time integration
The numerical parameters is given and the initial conditions are written in the output directory.

.. code-block:: python
  
  outf = 100                                                  # number of iterations between output files
  dt   = 1e-3                                                 # time step
  t    = 0                                                    # initial time
  tEnd = 1.5                                                  # final time
  ii   = 0                                                    # interator
  tic  = time.process_time()
  init_disk.write_mig(outputdir, t)

The computational loop can start.
The external forces given to the particles has to be given at each time step.

.. code-block:: python
  
  while t < tEnd:
      forces = g*np.pi*init_disk.r()**2*rhop
      time_integration.iterate(None, init_disk, dt, 2, 2e-6, forces)
      t  += dt
      ii += 1
      if ii%outf==0:
          init_disk.write_mig(outputdir, t)
      print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.process_time() - tic))
  t += dt

From disks to bodies
--------------------
We create a new problem for the bodies.
The option fixed_contact_geometry has to be put to zero to solve bodies dynamics.

.. code-block:: python
  
  p = scontact.ParticleProblem(2)
  p.set_friction_coefficient(friction)
  p.set_fixed_contact_geometry(0)

The open box is styill defined using three segments.

.. code-block:: python
  
  p.add_boundary_segment([0, 0], [0, h], "left")
  p.add_boundary_segment([w, h], [w, 0], "right")
  p.add_boundary_segment([w, 0], [0, 0], "bottom")

Each body of this new problem will have the position of the initial problem

.. code-block:: python
  
  p_bodies = np.asarray([], int)
  for i,body in enumerate(bodies):
      rpi      = init_disk.r()[i]/n_p_b
      x,y      = np.where(model)
      x        = 2*np.stack([x,y],axis=1)
      xb       = rpi*x
      radii    = np.repeat(rpi, n_p_b)
      masses   = radii**2*np.pi*rhop
      p_body   = p.add_particle_body(xb,radii,masses)
      p_bodies = np.append(p_bodies, p_body)
  p.body_position()[:] = init_disk.body_position()[:]
  p.body_theta()[p_bodies,:] = -angle

First, we will let the bodies falling without any rotation.
To do so, body inverse inertias are set to zero.

.. code-block:: python
  
  body_iinvertia = p.body_invert_inertia().copy()
  p.body_invert_inertia()[p_bodies] = 0
  

The numerical parameters are modified and the computational loop can start.

.. code-block:: python
  
  outf = 100
  dt   = 2.5e-4
  tEnd = 2
  while t < tEnd:
      forces = g*np.pi*p.r()**2*rhop
      time_integration.iterate(None, p, dt, 2, 2e-6, forces)
      t  += dt
      ii += 1
      if ii%outf==0:
          p.write_mig(outputdir, t)
      print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.process_time() - tic))
  t += dt
  

Then the rotation are activated by restoring the body inertias :

.. code-block:: python
  
  p.body_invert_inertia()[:] = body_iinvertia

The numerical parameters are modified and the computational loop can start.

.. code-block:: python
  
  outf = 100
  dt   = 2.5e-4
  tEnd = 3.5
  while t < tEnd:
      forces = g*np.pi*p.r()**2*rhop
      time_integration.iterate(None, p, dt, 2, 2e-6, forces)
      t  += dt
      ii += 1
      if ii%outf==0:
          p.write_mig(outputdir, t)
      print("ii : %4d"%ii)
  t += dt

The boundaries is then modified to describe the same domain as the fluid

.. code-block:: python
  
  w, h = width, height
  p.body_position()[:3] = 0                               # restore initial boundary body position
  p.segments()[0]["relative_position"]  = np.array([[0,0],[0,h]])         # left
  p.segments()[1]["relative_position"]  = np.array([[w,h],[w,0]])         # right
  p.segments()[2]["relative_position"]  = np.array([[w,0],[0,0]])         # bottom
  p.add_boundary_segment([0,h], [w,h], "top")

Mesh generation
----------------
The domain is created through the python GMSH api is loaded to generate the mesh and extract the boundaries.
The refinement is given by the function set_size_call_back. In this example, a constant mesh size is chosen.

.. code-block:: python
  
  mesh_size = 0.01                                # mesh size
  def gen_mesh(width, height, mesh_size, origin=np.array([0,0])):
      origin = np.asarray(origin)
      gmsh.model.add("box")
      gmsh.model.occ.add_rectangle(origin[0], origin[1], 0, width, height)
      gmsh.model.occ.synchronize()
      def get_line(x0, x1, eps =1e-6):
          r = gmsh.model.get_entities_in_bounding_box(x0[0]-eps, x0[1]-eps, -eps, x1[0]+eps, x1[1]+eps, eps, 1)
          return list(tag for dim, tag in r)
      h, w = height, width
      gmsh.model.add_physical_group(1,get_line(origin+[0,0], origin+[w,0]), name="Bottom")
      gmsh.model.add_physical_group(1,get_line(origin+[0,h], origin+[w,h]), name="Top")
      gmsh.model.add_physical_group(1,get_line(origin+[0,0], origin+[0,h]), name="Left")
      gmsh.model.add_physical_group(1,get_line(origin+[w,0], origin+[w,h]), name="Right")
      gmsh.model.add_physical_group(2,[1], name="domain")
      gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, lc : mesh_size)
      gmsh.model.mesh.generate(2)
  
  gen_mesh(w, h, mesh_size)
  

Fluid Problem
-------------
The fluid is described by its dimension, the external volume force applied, its dynamic viscosities and its densities.
The 2-fluid model is deduced from the length of the densities and viscosities provided.

.. code-block:: python
  
  rho0 = 1                                                # fluid density
  nu0  = 1e-5                                             # kinematic viscosity
  rho1 = 1000                                             # fluid density
  nu1  = 1e-6                                             # kinematic viscosity
  mu0  = rho0*nu0                                         # dynamic viscosity
  mu1  = rho1*nu1                                         # dynamic viscosity
  f   = fluid.FluidProblem(2,g,[mu0, mu1],[rho0, rho1], model_b=1)
  
  ## %%
  # The mesh can be loaded through the GMSH api if None is given or through a mesh filename.
  # The boundary conditions are described by their physical name.
  # To describle fully the pressure, its mean value over all the nodal values is impose to zero.
  f.load_msh(None)
  f.set_wall_boundary("Bottom")
  f.set_wall_boundary("Left")
  f.set_wall_boundary("Right")
  f.set_wall_boundary("Top")
  f.set_mean_pressure(0)

The initial air concentration is described by setting the concentration to 1 in the half upper domain.

.. code-block:: python
  
  x_cg     = f.coordinates()
  c0       = np.zeros_like(x_cg[:,0])
  x_up     = x_cg[:,1] > h/2
  c0[x_up] = 1
  f.set_concentration_cg(c0)
  

FEM-DEM coupling
----------------
The presence of particles is given to the fluid through the fluid volume fraction, i.e. the porosity and through a drag parametrization.
All the relevant informations needed for the fluid is given by :

.. code-block:: python
  
  f.set_particles(p.delassus(), p.volume(), p.position(), p.velocity(), p.contact_forces())
  

Time integration
----------------
The numerical parameters is given and the initial conditions are written in the output directory.

.. code-block:: python
  
  outf = 25                                               # number of iterations between output files
  dt   = 5e-4                                             # time step
  tEnd = 10                                               # final time
  ii   =  0                                               # initial iteration
  p.write_mig(outputdir, t)
  f.write_mig(outputdir, t)
  tic = time.process_time()
  

The computational loop can start.
The particles phase will use sub-iterations to keep a stable method. The minimal number of subiterations is given by min_nsub.
The external forces given to the particles has to be given at each time step.
To solve the two scales, first the fluid is solved then the particles are resolved.
A predictor-corrector method can also be used by imposing the flag to True.

.. code-block:: python
  
  while t < tEnd :
      forces = np.pi*p.r()**2*rhop * g
      time_integration.iterate(f, p, dt, min_nsub=2, external_particles_forces=forces, use_predictor_corrector=False)
      t  += dt
      ii += 1
      # Output files writting
      if ii%outf == 0 :
          p.write_mig(outputdir, t)
          f.write_mig(outputdir, t)
      print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.process_time() - tic)
