Download :download:`this example <../../examples/drop.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 fall of a cloud made of particles in a viscous 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 ctypes
  import numpy as np
  from scipy.spatial import KDTree
  from migflow import fluid
  from migflow import scontact
  from migflow import time_integration
  
  use_callback = True
  

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

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

Mesh generation
---------------
The mesh geometry is created with the python GMSH api.

.. code-block:: python
  
  height = 2                                      # domain height
  width = 0.1                                     # domain width
  origin = [-width/2, -height/2]                  # leftmost bottom corner
  gmsh.model
  def gen_geometry(width, height, 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])+get_line(origin+[w,0], origin+[w,h]), name="Lateral")
      gmsh.model.add_physical_group(2,[1], name="domain")

The mesh size can be chosen based on three approaches; based on a computed field, based on a callback or based on a uniform mesh size.

.. code-block:: python
  
  

The mesh size is computed based on the gradient variation. The field is then given to gmsh which will generate the mesh.

.. code-block:: python
  
  def gen_mesh_field(f):
      lcmin = 0.0015/4
      lcmax = 0.015/8   
      grad = f.fields_gradient()
      grad_u, grad_v, grad_p = grad[:,0,:], grad[:,1,:], grad[:,2,:]
      grad_v = np.linalg.norm(grad_v, axis=1)
      grad_v_min = np.min(grad_v)
      grad_v_max = np.max(grad_v)
      lv = (grad_v_max - grad_v_min)/(grad_v-grad_v_min) * lcmin
      size = np.maximum(lv, lcmax)
  
      x = f.coordinates()[f.elements()]
      size_view = gmsh.view.add("size")
      data = np.c_[x.swapaxes(1,2).reshape(-1,9), size[f.elements()]]
      size_field = gmsh.model.mesh.field.add("PostView")
      gmsh.model.mesh.field.setNumber(size_field, "ViewTag", size_view)
      gmsh.model.mesh.field.setAsBackgroundMesh(size_field)
      gmsh.view.add_list_data(size_view, "ST", data.shape[0], data.reshape(-1))
      gmsh.model.mesh.clear()
      gmsh.model.mesh.generate(2)
      gmsh.view.remove(size_view)
  

The mesh size is computed based on a callback function. Here a refinement is done close to the particles position.

.. code-block:: python
  
  def gen_mesh_callback(xp):
      tree = KDTree(xp)
      def size_f(x):
          dist, _ = tree.query(x[:,:2])
          distmin = 0.01
          lcmin = 0.0015
          lcmax = 0.015
          alpha = np.clip((dist[0]-distmin)/0.1, 0, 1)
          size = lcmin*(1-alpha) + lcmax*alpha
          return size
      gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, lc: size_f(np.array([[x,y]])))
      gmsh.model.mesh.clear()
      gmsh.model.mesh.generate(2)
      gmsh.model.mesh.remove_size_callback()
  

The mesh is generated based on a constant mesh size.

.. code-block:: python
  
  def gen_mesh_uniform(size):
      gmsh.model.mesh.clear()
      gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, lc: size)
      gmsh.model.mesh.generate(2)
      gmsh.model.mesh.remove_size_callback()
  

Particle Problem
----------------
The particle problem is created and its dimension is provided.

.. code-block:: python
  
  p = scontact.ParticleProblem(2)
  

The particle properties are defined.
All the particles are assumed to be spherical.

.. code-block:: python
  
  r    = 25e-6                                       # particles radius
  rhop = 2450                                         # particles density
  compacity = 0.2                                     # solid volume fraction in the drop
  rout = 3.3e-3                                       # drop radius
  

The particle positions are initialised randomly in a circular domain to obtain a given compacity.

.. code-block:: python
  
  def genInitialPosition(p, r, rout, rhop, compacity, origin) :
      """Set all the particles centre positions and create the particles objects
      
      Keyword arguments:    
      p -- Particle Problem
      r -- max radius of the particles
      rout -- outer radius of the cloud
      rhop -- particles density        
      compacity -- initial compacity in the cloud
      """
      # Space between the particles to obtain the expected compacity
      N = compacity*(rout/r)**2
      bodies = np.asarray([], int)
      while p.n_particles() < N:
          xyp = np.random.uniform(-rout, rout, 2)
          if np.hypot(xyp[0], xyp[1]) < rout:
              if p.n_particles() == 0:
                  d = 1
              else:
                  d = np.min(np.linalg.norm(p.position()-xyp[None,:],axis=1))
              if d > 2.1*r:
                  body = p.add_particle(xyp, r, r**2 * np.pi * rhop)
                  bodies = np.append(bodies,body)
      # Shift of the particles to the top of the box
      p.body_position()[bodies, :] += origin
      return bodies
  bodies = genInitialPosition(p, r, rout, rhop, compacity, np.array([0,1.9*height/4]))
  
  

Fluid Problem
-------------
The fluid is described by its dimension, the external volume force applied, its dynamic viscosity and its density.

.. code-block:: python
  
  g   = np.array([0,-9.81])                           # gravity
  rho = 1030                                          # fluid density
  nu  = 1.17/(2*rho)                                  # kinematic viscosity
  mu  = rho*nu                                        # dynamic viscosity
  f   = fluid.FluidProblem(2,g,[nu*rho],[rho], usolid=True)
  

The mesh is created with a refinment close to the particles position.
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
  
  gen_geometry(width, height, origin)
  if use_callback:
      gen_mesh_callback(p.position()[p.r()[:,0]>0])
  else:
      gen_mesh_uniform(0.015)
      f.load_msh(None)
  p.load_msh_boundaries(None, ["Bottom","Top","Lateral"])
  ## %%
  # 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 fully describe the pressure, its mean value over all the nodal values is set to zero.
  f.load_msh(None)
  f.set_wall_boundary("Bottom")
  f.set_wall_boundary("Top")
  f.set_wall_boundary("Lateral")
  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 are defined given and the initial conditions are written in the output directory.

.. code-block:: python
  
  outf = 10                                           # number of iterations between output files
  remeshing=20                                         # time step to adapt the mesh
  dt   = 5e-3                                         # time step
  tEnd = 50                                           # final time
  t    =  0                                           # initial time
  ii   =  0                                           # initial iteration
  

Write chosen fields into the output
-----------------------------------

.. code-block:: python
  
  def get_fields(fluid):
      y = fluid.coordinates_fields()[fluid.pressure_index()][:,1]
      y = y.reshape(-1,1)
      p1_element = fluid.get_p1_element()
      return {"pressure": (fluid.pressure(), p1_element),
              "velocity": (fluid.velocity(), p1_element),
              "porosity": (fluid.porosity(), p1_element),
              "dynamic_pressure":(fluid.pressure()-rho*g[1]*y, p1_element)}
  
  p.write_mig(outputdir, t)
  f.write_mig(outputdir, t, get_fields(f))

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 have to be given at each time step.
A predictor-corrector method can also be used by setting the flag to True.

.. code-block:: python
  
  tic = time.process_time()
  forces = g*np.pi*p.r()**2*rhop
  mass = np.pi*p.r()**2*rhop
  while t < tEnd :
      print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.process_time() - tic))
      f.set_particles(p.delassus(), p.volume(), p.position(), p.velocity(), p.omega(), p.contact_forces())
      tic = time.process_time()
      f.implicit_euler(dt, check_residual_norm=1)
      print("implicit : ", time.process_time()-tic)
      fext = f.compute_node_force(dt) + g*mass
      tic = time.process_time()
      time_integration._advance_particles(p, fext, dt, 1, 1e-3*r)
      print("NSCD : ", time.process_time()-tic)
  
      # time_integration.iterate(f, p, dt, min_nsub=2, external_particles_forces=forces, use_predictor_corrector=False)
      if ((ii+1)%remeshing==0):
          number_p = f.n_particles
          position_p = f.particle_position()
          volume_p = f.particle_volume()
      if (ii%remeshing==0 and ii != 0):
          if use_callback:
              gen_mesh_callback(p.position()[p.r()[:,0]>0])
          else:
              gen_mesh_field(f)
          f.adapt_mesh(old_n_particles=number_p, old_particle_position=position_p, old_particle_volume=volume_p)
          f.set_particles(p.delassus(), p.volume(), p.position(), p.velocity(), p.omega(), p.contact_forces())
      t += dt
      # Output files writting
      if ii%outf == 0 :
          p.write_mig(outputdir, t)
          f.write_mig(outputdir, t, get_fields(f))
      ii += 1
