Download :download:`this testcase <3d_drop_mono.py>`.

Three-dimensional fall of a single cloud made of particles in a viscous fluid
=============================================================================
This examples illustrates how to define the fluid problem and the particles one.
Setup is based on "Implementation of an unresolved stabilised FEM–DEM model to solve immersed granular flows" by M.Constant
We consider the Set 3 case.
All the relevant mesh information is directly loaded from the GMSH api.

Keywords
--------
FEM, DEM, Sedimentation, Mesh Adaptation

.. code-block:: python
  
  import gmsh
  import os, time, shutil, sys
  import numpy as np
  from scipy.spatial import KDTree
  from migflow import fluid, scontact
  
  np.random.seed(0)
  use_diagonal_offset = False
  

Output Directory
----------------
The output directory is (re)created

.. code-block:: python
  
  outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
  shutil.rmtree(outputdir, ignore_errors=True)
  os.makedirs(outputdir)
  

Physical Parameters
-------------------

.. code-block:: python
  
  rout = 2.7e-3  # drop radius
  r = 25e-6  # particles radius
  compacity = 0.03  # solid volume fraction in the drop
  g = np.array([0, -9.81, 0])  # gravity
  rho = 1220  # fluid density
  rhop = 2400  # particles density
  nu = 1e-4  # kinematic viscosity
  mu = rho * nu  # dynamic viscosity
  

Geometrical parameters
----------------------

.. code-block:: python
  
  height = 1  # domain height
  width = 10 * rout  # domain width
  depth = 10 * rout  # domain width
  origin = [-width / 2, -height / 2, -depth / 2]  # leftmost bottom corner
  
  

Mesh generation
---------------
The mesh is created with the python GMSH api.
The refinement is given by the function set_size_call_back. In this example, a refinement close to the the particles position is chosen.

.. code-block:: python
  
  def gen_geometry(width, height, origin=np.array([0, 0, 0])):
      origin = np.asarray(origin)
      gmsh.model.add("box")
      gmsh.model.occ.add_box(origin[0], origin[1], origin[2], width, height, depth)
      gmsh.model.occ.synchronize()
  
      def get_surface(x0, x1, eps=1e-6):
          r = gmsh.model.get_entities_in_bounding_box(
              x0[0] - eps,
              x0[1] - eps,
              x0[2] - eps,
              x1[0] + eps,
              x1[1] + eps,
              x1[2] + eps,
              2,
          )
          return list(tag for dim, tag in r)
  
      h, w, d = height, width, depth
      lateral = (
          get_surface(origin + [0, 0, 0], origin + [0, h, d])
          + get_surface(origin + [0, 0, 0], origin + [w, h, 0])
          + get_surface(origin + [w, 0, 0], origin + [w, h, d])
          + get_surface(origin + [0, 0, d], origin + [w, h, d])
      )
      gmsh.model.add_physical_group(
          2, get_surface(origin + [0, h, 0], origin + [w, h, d]), name="Top"
      )
      gmsh.model.add_physical_group(
          2, get_surface(origin + [0, 0, 0], origin + [w, 0, d]), name="Bottom"
      )
      gmsh.model.add_physical_group(2, lateral, name="Lateral")
      gmsh.model.add_physical_group(3, [1], name="Domain")
  
  
  def gen_mesh_callback(xp):
      tree = KDTree(xp)
  
      def size_f(x):
          dist, _ = tree.query(x[:, :3])
          distmin = 0.005
          mylc = 0.002
          alpha = np.clip((dist[0] - distmin) / 0.1, 0, 1)
          size = mylc * (1 - alpha) + 10 * mylc * alpha
          return size
  
      gmsh.model.mesh.set_size_callback(
          lambda dim, tag, x, y, z, lc: size_f(np.array([[x, y, z]]))
      )
      gmsh.model.mesh.clear()
      gmsh.model.mesh.generate(3)
      gmsh.model.mesh.remove_size_callback()
  
  

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

.. code-block:: python
  
  def genInitialPosition(p: scontact.ParticleProblem, 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
      """
      N_p = int(compacity * (rout / r) ** 3)
      bodies = np.asarray([], int)
      n = 0
      while n < N_p:
          print("...add particle : ", n, "/", N_p, end="\r")
          xyp = np.random.uniform(-rout, rout, 3)
          if np.linalg.norm(xyp) < 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, 4 / 3 * r**3 * np.pi * rhop)
                  bodies = np.append(bodies, body)
                  n += 1
      # Shift of the particles to the top of the box
      actual_compacity = np.sum(p.volume()) / ((4 / 3) * np.pi * rout**3)
      print(f"Actual compacity: {actual_compacity:.3f} ({p.n_particles()} particles)")
      p.body_position()[bodies, :] += origin
      return bodies
  
  

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

.. code-block:: python
  
  p = scontact.ParticleProblem(3)
  offset0 = np.array([-2 * rout if use_diagonal_offset else 0, 1.9 * height / 4, 0])
  genInitialPosition(p, r, rout, rhop, compacity, offset0)
  

Mesh generation
---------------
The mesh is created with the python GMSH api.
The refinement is given by the function set_size_call_back. In this example, a refinement close to the the particles position is chosen.

.. code-block:: python
  
  gen_geometry(width, height, origin)
  gen_mesh_callback(p.position()[p.r()[:, 0] > 0])
  # The mesh can be loaded through the GMSH api if None is given or through a mesh filename.
  p.load_msh_boundaries(None, ["Bottom", "Top", "Lateral"])
  

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

.. code-block:: python
  
  f = fluid.FluidProblem(3, g, nu * rho, rho)
  f.load_msh(None)
  f.set_wall_boundary("Bottom")
  f.set_wall_boundary("Top")
  f.set_wall_boundary("Lateral")
  f.set_mean_pressure(0)
  
  

Simulation Parameters
---------------------

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

Write specific fields into the output
-------------------------------------

.. code-block:: python
  
  def dynamic_pressure(f, rho, g):
      y = f.coordinates_fields()[f.pressure_index()][:, 1]
      return f.pressure() - rho * g[1] * y
  
  
  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),
          "u_solid": (fluid.u_solid(), p1_element),
          "dynamic_pressure": (fluid.pressure() - rho * g[1] * y, p1_element),
      }
  
  

Computational Loop
------------------

.. code-block:: python
  
  tic = time.process_time()
  while t < tEnd:
      print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.process_time() - tic))
      if ii % remeshing == 0 and ii != 0:
          gen_mesh_callback(p.position()[p.r()[:, 0] > 0])
          f.adapt_mesh()
      f.set_particles(
          p.delassus(),
          p.volume(),
          p.position(),
          p.velocity(),
          p.omega(),
          p.contact_forces(),
      )
      if ii % outf == 0:
          f.write_mig(outputdir, t, get_fields(f))
          p.write_mig(outputdir, t)
      f.implicit_euler(dt)
      forces = 4 / 3 * g * np.pi * p.r() ** 3 * rhop + f.compute_node_force()
      p.iterate(dt, forces, tol=1e-3 * r, force_motion=1)
      t += dt
      ii += 1
