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

Granular Column Collapse in Fluid
=================================
This test case simulates the collapse of a granular column in a viscous fluid.
The column of grains can either be imported from a pre-deposited configuration
(see ``../depot/depot.py``) or generated automatically using a regular
hexagonal packing. The case demonstrates the coupling between a dense
granular phase and a single-phase incompressible fluid.

Keywords
--------
FEM, DEM, Unresolved, Friction

Description
-----------
The granular column initially rests in hydrostatic equilibrium inside a
rectangular box filled with water. Upon release, the grains collapse and
form a dense avalanche driven by gravity and resisted by viscous drag.


.. code-block:: python
  
  import os, sys, shutil, subprocess
  import time
  import numpy as np
  from migflow import fluid, scontact, time_integration
  
  np.random.seed(0)
  

Output Directory
----------------
Create an output directory for results and generate a 2D mesh using Gmsh.

.. code-block:: python
  
  outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
  geomesh_filename = f"mesh_unresolved.geo"
  mesh_filename = f"{outputdir}/mesh.msh"
  shutil.rmtree(outputdir, ignore_errors=True)
  os.makedirs(outputdir)
  subprocess.call(["gmsh", "-2", geomesh_filename, "-o", mesh_filename])
  

Hexagonal Packing Function
--------------------------
Define a helper function to create a regular hexagonal packing of grains when
no pre-deposited configuration is available.

.. code-block:: python
  
  use_pre_deposit = False
  
  
  def hexagonal_packing(p, x0, y0, lx, ly, rmin, rmax):
      """Create a hexagonal particle arrangement inside a rectangular region.
  
      Parameters
      ----------
      p : migflow.scontact.ParticleProblem
          Particle problem instance.
      x0, y0 : float
          Bottom-left corner of the region.
      lx, ly : float
          Width and height of the region.
      rmin, rmax : float
          Minimum and maximum particle radii.
      """
      for x in np.arange(x0 + rmax, x0 + lx - rmax, 2 * rmax):
          for y in np.arange(y0 + rmax, y0 + ly - rmax, 2 * (3**0.5) * rmax):
              z = np.random.random((1))
              r = rmin + (rmax - rmin) * z
              p.add_particle((x, y), r, np.pi * r**2 * rhop)
      for x in np.arange(2 * rmax + x0, x0 + lx - rmax, 2 * rmax):
          for y in np.arange(
              y0 + (3**0.5) * rmax + rmax, y0 + ly - rmax, 2 * (3**0.5) * rmax
          ):
              z = np.random.random((1))
              r = rmin + (rmax - rmin) * z
              p.add_particle((x, y), r, np.pi * r**2 * rhop)
  
  

Physical Parameters
-------------------
The simulation parameters are defined below, including fluid viscosity,
densities, gravity, and the typical particle radius.

.. code-block:: python
  
  g = np.array([0, -9.81])  # gravity [m/s²]
  rhop = 1500  # particle density [kg/m³]
  rho = 1000  # fluid density [kg/m³]
  nu = 1e-6  # fluid kinematic viscosity [m²/s]
  mu = rho * nu  # fluid dynamic viscosity [Pa·s]
  r = 2e-3  # mean particle radius [m]
  h = 0.2  # domain height [m]
  lx = 0.1  # granular column width [m]
  ly = 0.165  # granular column height [m]
  
  

Particle Problem
----------------
The particle problem is created and the boundary geometry is loaded.
Depending on ``use_pre_deposit``, the initial configuration is either
imported from a previous *depot* testcase or generated on the fly.

.. code-block:: python
  
  p = scontact.ParticleProblem(2)
  p.set_fixed_contact_geometry(0)
  p.load_msh_boundaries(mesh_filename, ["Top", "Bottom", "Left", "Right"])
  
  if use_pre_deposit:
      # Load compact column from depot testcase
      p1 = scontact.ParticleProblem(2)
      p1.read_mig("../depot/output", 10)
      p.add_particles(p1.position(), p1.r(), p1.mass(), v=p1.velocity())
  else:
      # Generate hexagonal packing in the left side of the domain
      eps = 1e-4 * r
      hexagonal_packing(p, eps, eps, lx, ly, 0.9 * r, 1.1 * r)
  
  # Set particle–particle and particle–wall friction
  p.set_friction_coefficient(0.5)
  p.write_mig(outputdir, 0)
  

Fluid Problem
-------------
The fluid problem is initialized using the same geometry and coupled to the
particle phase. All walls are considered no-slip, and the mean pressure is
set to zero.

.. code-block:: python
  
  f = fluid.FluidProblem(2, g, mu, rho)
  f.load_msh(mesh_filename)
  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)
  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
  
  i = 0
  outf = 10  # number of iterations between outputs
  dt = 1e-3  # time step [s]
  t = 0  # initial time [s]
  tEnd = 5  # final time [s]
  mass = np.pi * p.r() ** 2 * rhop
  tic = time.time()
  while t < tEnd:
      print(f"{i = :4d} ----- {t = :g}")
      if i % outf == 0:
          p.write_mig(outputdir, t)
          f.write_mig(outputdir, t)
      time_integration.iterate(
          f, p, dt, min_nsub=1, external_particles_forces=g * mass, contact_tol=1e-4 * r
      )
      t += dt
      i += 1
  

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

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

