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

2D Granular column collapse onto a fluid bed
============================================
This example simulates the collapse of a granular column onto a fluid bed using the PFEM-DEM approach.
It showcases the interaction between granular materials and fluids, capturing the dynamics of the collapse and subsequent fluid-grain interactions.
The aim is to visualise the type of free surface wave generated by the impact of the grains on the fluid bed.
This test case follows the setup described in:
- W. Sarlin, C. Morize, A. Sauret, P. Gondret, Nonlinear regimes of tsunami waves generated by a granular collapse, J. Fluid. Mech. 919 (2021) R6.

.. code-block:: python
  
  
  # Keywords
  # --------
  # PFEM, fluid-grains, granular column collapse, 2D
  #
  # Description
  # -----------
  # The test demonstrates:
  # - How to perform a PFEM-DEM simulation with moving boundaries.
  import numpy as np
  import sys
  import shutil
  import os
  

Gmsh pfem branch
----------------
This test requires the alphashapes branch of GMSH to use the Particle Finite Element (PFEM) module.
PFEM is a lagrangian method where the mesh is regenerated at each time step, making it suitable for problems with large deformations and free surfaces.

.. code-block:: python
  
  sys.path.insert(0, os.environ["GMSH_PFEM_DIR"])
  
  import gmsh
  from migflow import fluid, time_integration, scontact, pfem
  

Output Directory
----------------
Create a clean output directory for simulation results.

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

Geometrical parameters and mesh generation
------------------------------------------

.. code-block:: python
  
  L = 2.0  # length of the domain
  H = 1.0  # height of the domain
  H_0_array = np.array([0.39, 0.29, 0.29])  # initial column heights
  L_0_array = np.array([0.145, 0.1, 0.05])  # initial column lengths
  h_0_array = np.array([0.04, 0.08, 0.25])  # initial fluid heights
  
  H_0 = H_0_array[0]
  L_0 = L_0_array[0]
  l_0 = L - L_0
  h_0 = h_0_array[0]
  t_gate = 0.01 * L_0  # gate thickness
  y_gate = h_0 + 1e-4  # initial gate position
  h_gate = y_gate + 1.2 * H_0  # gate height
  u_gate = 1  # m/s; gate upward velocity
  

Mesh parameters
---------------

.. code-block:: python
  
  geo_mesh_size = L / 10
  mesh_size = 0.01 * l_0
  smax = mesh_size
  smin = 0.2 * mesh_size
  dmax = 6 * smax
  alpha = 1.3
  

Initial fluid mesh
------------------

.. code-block:: python
  
  gmsh.model.add("ModelInit")
  gmsh.model.occ.addRectangle(L_0, 0, 0, l_0, h_0)
  gmsh.model.occ.synchronize()
  gmsh.model.mesh.setSizeCallback(lambda *args: mesh_size)
  gmsh.model.mesh.generate(2)
  nodeTags, coords, _ = gmsh.model.mesh.getNodes()
  

Solid domain
------------

.. code-block:: python
  
  gmsh.model.add("ModelGeo")
  
  
  def set_geo_mesh(y_gate=1e-4):
      name = gmsh.model.getCurrent()
      gmsh.model.setCurrent("ModelGeo")
      dimtag = gmsh.model.getEntities(2)
      gmsh.model.occ.remove(dimtag, True)
      box = gmsh.model.occ.addRectangle(0, 0, 0, L, H)
      step = gmsh.model.occ.addRectangle(0, 0, 0, L_0, h_0)
      gate = gmsh.model.occ.addRectangle(L_0, y_gate, 0, t_gate, h_gate)  # gate
      gmsh.model.occ.cut([(2, box)], [(2, gate), (2, step)])
      gmsh.model.occ.synchronize()
      gmsh.model.mesh.setSizeCallback(lambda *args: geo_mesh_size)
      gmsh.model.mesh.generate(2)
      geoEntities = gmsh.model.getEntities(1)
      gmsh.model.setCurrent(name)
      return geoEntities
  
  
  geoEntities = set_geo_mesh(y_gate)
  gmsh.model.addPhysicalGroup(1, [16], -1, "bottom")
  gmsh.model.addPhysicalGroup(1, [15], -1, "right")
  gmsh.model.addPhysicalGroup(1, [13], -1, "left")
  gmsh.model.addPhysicalGroup(1, [7], -1, "step_top")
  gmsh.model.addPhysicalGroup(1, [6], -1, "step_right")
  gmsh.model.addPhysicalGroup(1, [9, 10, 11, 12], -1, "gate")
  

PFEM mesh and size fields
-------------------------

.. code-block:: python
  
  gmsh.model.add("ModelFluid")
  alphaDomainTag = gmsh.model.addDiscreteEntity(2, -1, [])
  for dim, tag in geoEntities:
      gmsh.model.addDiscreteEntity(dim, tag, [])
  alphaBoundaryTag = gmsh.model.addDiscreteEntity(1, -1, [])
  gmsh.model.mesh.addNodes(2, alphaDomainTag, nodeTags, coords)
  gmsh.model.addPhysicalGroup(1, [16], -1, "bottom")
  gmsh.model.addPhysicalGroup(1, [15], -1, "right")
  gmsh.model.addPhysicalGroup(1, [13], -1, "left")
  gmsh.model.addPhysicalGroup(1, [7], -1, "step_top")
  gmsh.model.addPhysicalGroup(1, [6], -1, "step_right")
  gmsh.model.addPhysicalGroup(1, [9, 10, 11, 12], -1, "gate")
  gmsh.model.addPhysicalGroup(1, [alphaBoundaryTag], -1, "freeSurface")
  gmsh.model.addPhysicalGroup(2, [alphaDomainTag], -1, "domain")
  

Size fields
-----------

.. code-block:: python
  
  sizeFieldConstant = gmsh.model.mesh.field.add("Box")
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "VIn", mesh_size)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "VOut", 10 * mesh_size)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "XMin", 0.0)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "XMax", L)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMin", 0.0)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMax", H)
  gmsh.model.mesh.field.setNumber(sizeFieldConstant, "Thickness", 0.001)
  
  sizeFieldDistFS = gmsh.model.mesh.field.add("AlphaShapeDistance")
  gmsh.model.mesh.field.setNumber(sizeFieldDistFS, "Tag", alphaBoundaryTag)
  gmsh.model.mesh.field.setNumber(sizeFieldDistFS, "SamplingLength", 0.25 * smin)
  sizeFieldRefine1 = gmsh.model.mesh.field.add("Threshold")
  gmsh.model.mesh.field.setNumber(sizeFieldRefine1, "InField", sizeFieldDistFS)
  gmsh.model.mesh.field.setNumber(sizeFieldRefine1, "SizeMin", smin)
  gmsh.model.mesh.field.setNumber(sizeFieldRefine1, "SizeMax", smax)
  gmsh.model.mesh.field.setNumber(sizeFieldRefine1, "DistMin", 0.0)
  gmsh.model.mesh.field.setNumber(sizeFieldRefine1, "DistMax", dmax)
  
  bndEntities = []
  sizeFieldBnd = []
  for b in bndEntities:
      sf = gmsh.model.mesh.field.add("AlphaShapeDistance")
      gmsh.model.mesh.field.setNumber(sf, "Tag", b)
      gmsh.model.mesh.field.setNumber(sf, "SamplingLength", 0.25 * smin)
      sizeFieldBnd.append(gmsh.model.mesh.field.add("Threshold"))
      gmsh.model.mesh.field.setNumber(sizeFieldBnd[-1], "InField", sf)
      gmsh.model.mesh.field.setNumber(sizeFieldBnd[-1], "SizeMin", smin)
      gmsh.model.mesh.field.setNumber(sizeFieldBnd[-1], "SizeMax", smax)
      gmsh.model.mesh.field.setNumber(sizeFieldBnd[-1], "DistMin", 0.0)
      gmsh.model.mesh.field.setNumber(sizeFieldBnd[-1], "DistMax", dmax)
  sizeFieldRefine = gmsh.model.mesh.field.add("Min")
  gmsh.model.mesh.field.setNumbers(
      sizeFieldRefine, "FieldsList", [sizeFieldRefine1] + sizeFieldBnd
  )
  gmsh.model.mesh.computeAlphaShape(
      2,
      alphaDomainTag,
      alphaBoundaryTag,
      "ModelGeo",
      alpha,
      sizeFieldConstant,
      sizeFieldRefine,
      False,
  )
  

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

.. code-block:: python
  
  g = np.array([0.0, -9.81])
  rho = 1000
  mu = 1e-3
  rho_bubble = 1
  sigma = 0.00728
  use_bubble = False
  

DEM parameters
--------------

.. code-block:: python
  
  volume_corr = 0.8  # volume correction on the grains to account for 2D-3D discrepancies
  rhop = 2500
  friction = 0.9  # 0.2
  friction_wall = 0.9
  friction_gate = 0.0
  r_particles = 0.0025
  n_particles = int(L_0 / (2 * r_particles) * H_0 / (2 * r_particles) * 1.5)
  

Time parameters
---------------

.. code-block:: python
  
  cfl = 0.1
  U = 10
  U_init = 1.0
  t = 0
  dt = smin / U * cfl
  settle_time = 0.03
  tEnd = 1.0 + settle_time
  outf = 10
  print(f"dt = {dt}, U = {U}")
  

Fluid problem
-------------

.. code-block:: python
  
  f = fluid.FluidProblem(2, g, mu, rho, advection=False)
  f.set_wall_boundary("bottom")
  f.set_wall_boundary("right")
  
  
  def gate_velocity(t):
      if t > settle_time:
          return u_gate
      return 0
  
  
  f.set_wall_boundary("gate", velocity=[0, gate_velocity(t)])
  f.set_wall_boundary("left")
  f.set_wall_boundary("step_top")
  f.set_wall_boundary("step_bottom")
  f.set_strong_boundary("bottom", velocity=[0, 0])
  f.set_strong_boundary("right", velocity_x=0)
  f.set_strong_boundary("left", velocity_x=0)
  f.set_strong_boundary("step_top", velocity_y=0)
  f.set_strong_boundary("step_right", velocity_x=0)
  f.set_open_boundary("freeSurface", pressure=0, viscous_flag=False)
  

DEM Problem
-----------

.. code-block:: python
  
  dem = scontact.ParticleProblem(2)
  outer_body = dem.add_body((0, 0), 0, 0)
  x0 = np.array([L_0, 0])
  x1 = np.array([L, 0])
  x2 = np.array([L, H])
  x3 = np.array([0, H])
  x4 = np.array([0, h_0])
  x5 = np.array([L_0, h_0])
  dem.add_segment_to_body(x0, x1, outer_body, "Wall", material="Wall")
  dem.add_segment_to_body(x1, x2, outer_body, "Wall", material="Wall")
  dem.add_segment_to_body(x2, x3, outer_body, "Wall", material="Wall")
  dem.add_segment_to_body(x3, x4, outer_body, "Wall", material="Wall")
  dem.add_segment_to_body(x4, x5, outer_body, "Wall", material="Wall")
  dem.add_segment_to_body(x5, x0, outer_body, "Wall", material="Wall")
  dem.add_particle_to_body(x0, 0, outer_body)
  dem.add_particle_to_body(x1, 0, outer_body)
  dem.add_particle_to_body(x2, 0, outer_body)
  dem.add_particle_to_body(x3, 0, outer_body)
  dem.add_particle_to_body(x4, 0, outer_body)
  dem.add_particle_to_body(x5, 0, outer_body)
  
  gate_body = dem.add_body((0, 0), 0, 0)
  x0 = np.array([L_0, y_gate])
  x1 = np.array([L_0 + t_gate, y_gate])
  x2 = np.array([L_0 + t_gate, h_gate])
  x3 = np.array([L_0, h_gate])
  dem.add_segment_to_body(x0, x1, gate_body, "Gate", material="Gate")
  dem.add_segment_to_body(x1, x2, gate_body, "Gate", material="Gate")
  dem.add_segment_to_body(x2, x3, gate_body, "Gate", material="Gate")
  dem.add_segment_to_body(x3, x0, gate_body, "Gate", material="Gate")
  dem.add_particle_to_body(x0, 0, gate_body)
  dem.add_particle_to_body(x1, 0, gate_body)
  dem.add_particle_to_body(x2, 0, gate_body)
  dem.add_particle_to_body(x3, 0, gate_body)
  dem.body_velocity()[gate_body, 1] = u_gate
  
  dem.set_fixed_contact_geometry(0)
  dem.set_friction_coefficient(10)
  
  dem.set_friction_coefficient(friction, "Grain", "Grain")
  dem.set_friction_coefficient(friction_gate, "Grain", "Gate")
  dem.set_friction_coefficient(friction_wall, "Grain", "Wall")
  

initial grains deposition
-------------------------

.. code-block:: python
  
  ratio = 1.2
  diameter_max = 2.0 * r_particles
  diameter_min = diameter_max / ratio
  np.random.seed(0)
  u = np.random.power(1, n_particles)  # probability density function [0,1]
  d = (
      diameter_max * diameter_min / (diameter_max + u * (diameter_min - diameter_max))
  )  # random particle diameter following the given "u" law
  radii = d / 2
  boundaries = np.array([[0, H], [0, h_0], [L_0, h_0], [L_0, H]])
  xp, radii = scontact.fastdeposit_2d(boundaries, radii, "deposit.svg")
  eps = 1e-4 * diameter_max
  for x, r in zip(xp, radii):
      if x[1] < h_0 + H_0:
          dem.add_particle(x + eps, r, np.pi * r**2 * rhop, "Grain")
  

Simulation Loop
---------------
Time integration of coupled fluid–particle motion.

.. code-block:: python
  
  i = 0
  gmsh.option.setNumber("General.Verbosity", 0)
  while t < tEnd:
      print(f"{i:4d}, {t:.6g}/{tEnd:.6g}, {dt:.6g}")
      nodetag, oelemtag, oparamcoord, _ = gmsh.model.mesh.computeAlphaShape(
          2,
          alphaDomainTag,
          alphaBoundaryTag,
          "ModelGeo",
          alpha,
          sizeFieldRefine,
          sizeFieldRefine,
          boundaryTolerance=0.001 * smin,
          usePreviousMesh=True,
      )
      oparamcoord = oparamcoord.reshape((-1, 3))
  
      ordered_node_tags = pfem.prepareMeshForMigflow(
          i, alphaDomainTag, f, nodetag, oelemtag, oparamcoord
      )
      f.set_particles(
          dem.delassus() * volume_corr,
          dem.volume() * volume_corr,
          dem.position(),
          dem.velocity(),
          dem.omega() * 0,
          dem.contact_forces(),
      )
  
      pfem.applySurfaceTension(f, sigma, False, g, rho_bubble)
  
      if i % outf == 0:
          f.write_mig(outputdir, t)
          dem.write_mig(outputdir, t)
  
      f.implicit_euler(dt)
  
      if t < settle_time:
          dem.body_velocity()[gate_body, 1] = 0
      else:
          dem.body_velocity()[gate_body, 1] = u_gate
  
      drag = f.compute_node_force()
      external_forces = drag + g * dem.r() ** 2 * np.pi * rhop * volume_corr
      time_integration._advance_particles(dem, external_forces, dt, 1, 1e-6)
  
      if y_gate > H / 3:
          dem.body_velocity()[gate_body, 1] = 0
      y_gate += dem.body_velocity()[gate_body, 1] * dt
      geoEntities = set_geo_mesh(y_gate)
  
      dx = np.zeros((f.coordinates().shape[0], 3))
      dx[:, :2] = f.velocity() * dt
      dx = dx.reshape((-1,))
      gmsh.model.mesh.advectMeshNodes(
          2,
          alphaDomainTag,
          alphaBoundaryTag,
          "ModelGeo",
          ordered_node_tags,
          dx,
          0.01 * smin,
      )
  
      _, newCoords, _ = gmsh.model.mesh.getNodes(2, alphaDomainTag)
  
      f.set_coordinates(newCoords)
  
      t += dt
      i += 1
  
      max_U = np.max(
          [
              U_init,
              np.max(np.linalg.norm(f.velocity(), axis=1)),
              np.max(np.linalg.norm(dem.velocity(), axis=1)),
          ]
      )
      dt = smin / max_U * cfl
  

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

 python3 -m migflow.plot.migplot output --actors fluid particles --fluid-field velocity --show-edges 1

