Download :download:`this example <../../examples/depot.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 particles sedimentation in fluid
==============================================
This examples illustrates how to define the fluid problem as well as the particles one.
All the relevant mesh information is directly loaded from the GMSH api.

.. code-block:: python
  
  
  import os
  import time
  import gmsh
  import shutil
  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-depot"
  shutil.rmtree(outputdir,ignore_errors=True)
  if not os.path.isdir(outputdir) :
      os.makedirs(outputdir)
  

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
  
  
  height = 0.6                                    # domain height
  width = 0.4                                     # domain width
  mesh_size = 0.02                                # mesh size
  origin = [-width/2, -height/2]                  # leftmost bottom corner
  
  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(width, height, mesh_size, [-width/2, -height/2])
  

Particle Problem
----------------
The particle problem based is created, its dimension has to be provided.
The boundaries are loaded either from a .msh file or from the current model loaded with the GMSH api (if None as a filename is given). 

.. code-block:: python
  
  p = scontact.ParticleProblem(2)
  p.load_msh_boundaries(None, ["Top", "Left","Right","Bottom"])
  

Let's define the particle properties and initialised them on a compact rectangular grid at the top of the domain.
All the particles are assumed to be spherical.

.. code-block:: python
  
  r    = 2.5e-3                                       # particles radius
  rhop = 1500                                         # particles density
  h    = 1.5e-1                                       # particles area height
  w    = 4e-1                                         # particles area widht
  eps  = 1e-8                                         # tolerance
  x    = np.arange(r-eps, w-r+eps, 2*r) + origin[0]
  y    = np.arange(-r, -h, -2*r) - origin[1]
  x, y = np.meshgrid(x,y)
  for xi,yi in zip(x.flat,y.flat):
      p.add_particle((xi, yi), r, r**2*np.pi*rhop)
  

Fluid Problem
-------------
The fluid is described by its dimension, the external volume force applied, its dynamic viscosity and its density.
Moreover, in order to have a stable method for compact configuation, an extra-diffusion is activated by the drag_in_stab option.

.. code-block:: python
  
  g   = np.array([0,-9.81])                           # gravity
  rho = 1000                                          # fluid density
  nu  = 1e-6                                          # kinematic viscosity
  mu  = rho*nu                                        # dynamic viscosity
  f   = fluid.FluidProblem(2,g,nu*rho,rho)
  

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.

.. code-block:: python
  
  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)
  

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.omega(), 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   = 1e-3                                         # time step
  tEnd = 10                                           # final time
  t    =  0                                           # intial 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
      # Output files writting
      if ii%outf == 0 :
          p.write_mig(outputdir, t)
          f.write_mig(outputdir, t)
      ii += 1
      print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.process_time() - tic))
