Download :download:`this testcase <2d_avalanche_semi_resolved.py>`.

Semi-Resolved Avalanche
=======================
This example simulates a two-dimensional dense granular avalanche using a
semi-resolved fluid–particle coupling. In this formulation, the mesh size
of the fluid is of the same order as the particle diameter, allowing the
hydrodynamic forces to be partially resolved around individual particles.

A set of polydisperse particles is first deposited inside a rectangular
domain. The particles interact through frictional contacts,
while a viscous incompressible fluid surrounds them and exerts
drag and buoyancy. The simulation illustrates the onset of an avalanche and
the interplay between solid contact chains and fluid stresses.

Keywords
--------
FEM, DEM, Semi-Resolved, Adaptation, Friction, Polydispersity

.. code-block:: python
  
  import sys, os, shutil
  import gmsh
  import numpy as np
  from scipy.spatial import KDTree
  

This testcase relies on the LMGC90 library

.. code-block:: python
  
  sys.path.insert(0, os.environ["LMGC90_DIR"])
  
  from pylmgc90 import pre
  from migflow import fluid, scontact, time_integration
  
  np.random.seed(0)
  
  

Output Directory
----------------
The output directory is created (deleted first if it already exists).

.. 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
-------------------
Fluid and particle properties are defined here. The particles are
denser than the fluid to initiate motion under gravity.

.. code-block:: python
  
  rho = 1000.0  # fluid density [kg/m³]
  mu = 1e-3  # fluid viscosity [Pa·s]
  gravity = np.array([0.0, -9.81])  # gravity vector [m/s²]
  rhop = 2500.0  # particle density [kg/m³]
  dmax = 1e-2  # maximum particle diameter [m]
  dmin = 5e-3  # minimum particle diameter [m]
  polydispersity = dmax / dmin  # size ratio
  

Geometrical Parameters
----------------------
The computational domain is a rectangular box.
The mesh is generated so that the element size is of the same order
as the particle diameter (semi-resolved scale).

.. code-block:: python
  
  height = 0.4  # domain height [m]
  width = 0.8  # domain width [m]
  mesh_size = 10 * dmin  # mesh size [m]
  origin = np.array([-width / 2, -height / 2])  # mesh origin [m]
  eps = 1e-8  # numerical tolerance [m]
  h0 = 0.1875  # initial particle bed height [m]
  w0 = 1e-1  # initial particle bed width [m]
  
  

Mesh Generation
---------------
The Gmsh model is generated with physical groups for the walls and domain.
These will be used later for setting boundary conditions in the fluid and
particle solvers.
The local mesh size is controlled using the distance to the particle cloud.
This ensures refinement near the granular region and coarser elements
elsewhere.

.. code-block:: python
  
  lc = 2 * dmax
  lcmin = lc / 2
  lcmax = lc
  distmin = 10 * dmin
  
  
  class MeshSize:
      def __call__(self, dim, tag, x, y, z, lc):
          dist, _ = self.tree.query([[x, y]])
          alpha = np.clip((dist[0] - distmin) / (1 * distmin), 0, 1)
          return lcmin * (1 - alpha) + lcmax * alpha
  
      def set_points(self, points):
          self.tree = KDTree(points)
  
  
  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(mesh_size)
      gmsh.model.mesh.generate(2)
  
  

Particle Problem
----------------
The particle assembly is initialized by randomly depositing polydisperse
grains inside a rectangular region. Each grain is represented by a circular
particle with its own radius and mass. Particles are deposited inside the box using the LMGC preprocessor.
This generates a dense initial packing.

.. code-block:: python
  
  p = scontact.ParticleProblem(2)
  n_particles_max = int(2e5)  # upper bound for number of particles
  u = np.random.power(1, n_particles_max)  # probability density function [0,1]
  d = (
      dmax * dmin / (dmax + u * (dmin - dmax))
  )  # random particle diameter following the given "u" law
  rp = d / 2
  [nb_deposited, x, rp] = pre.depositInBox2D(rp, w0, h0)
  xp = x.reshape(-1, 2)
  xp += origin[None, :]
  for xi, ri in zip(xp, rp):
      p.add_particle(xi, ri, np.pi * ri**2 * rhop, "particle")
  

Mesh Adaptation
---------------
The adaptive mesh size field is built from the particle positions and used
to generate the mesh.

.. code-block:: python
  
  mesh_size = MeshSize()
  mesh_size.set_points(p.position()[p.r()[:, 0] > 0])
  gen_mesh(width, height, mesh_size, origin)
  

Boundary and Friction Conditions
--------------------------------
Wall boundaries are imported from the mesh and assigned frictional contact
properties. Particle–wall and particle–particle friction coefficients are
also defined.

.. code-block:: python
  
  p.load_msh_boundaries(None, ["Top", "Left", "Right", "Bottom"], material="wall")
  p.set_friction_coefficient(1.3, "wall", "particle")
  p.set_friction_coefficient(0.4, "particle", "particle")
  
  

Fluid Problem
-------------
The fluid problem is created and initialized from the Gmsh mesh. No-slip
boundary conditions are imposed on all walls. The particle information
(positions, velocities, contact forces) is passed to the fluid solver for
coupling.

.. code-block:: python
  
  f = fluid.FluidProblem(2, gravity, mu, rho)
  f.load_msh(None)
  f.set_wall_boundary("Bottom", velocity=[0, 0])
  f.set_wall_boundary("Left", velocity=[0, 0])
  f.set_wall_boundary("Right", velocity=[0, 0])
  f.set_wall_boundary("Top", velocity=[0, 0])
  f.set_mean_pressure(0)
  f.set_particles(
      p.delassus(), p.volume(), p.position(), p.velocity(), p.omega(), p.contact_forces()
  )
  

Time Integration
----------------
The simulation runs for a short time interval with a fixed time step.
At each iteration, the fluid phase is solved and the fluid-grain forces are applied to the particles,
and then system evolves according to the contact dynamics solver.
The simulation runs until the avalanche stabilizes.

.. code-block:: python
  
  t = 0
  i = 0
  dt = 1e-3
  tEnd = 4
  outf = 50
  mass = np.pi * p.r() ** 2 * rhop
  while t < tEnd:
      print(f"{i = :4d} ----- {t = :g}")
      if i % outf == 0:
          f.write_mig(outputdir, t)
          p.write_mig(outputdir, t)
      f.implicit_euler(dt, check_residual_norm=1)
      fext = gravity * mass + f.compute_node_force()
      time_integration._advance_particles(p, fext, dt, 1, 1e-3 * dmin)
      f.set_particles(
          p.delassus(),
          p.volume(),
          p.position(),
          p.velocity(),
          p.omega(),
          p.contact_forces(),
      )
      t += dt
      i += 1
  

Plot
----
.. code-block:: shell

 python3 -m migflow.plot.migplot output --actors fluid particles --fluid-field pressure

