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

Bidimensional Particle Sedimentation (Grid Packing)
===================================================
This example illustrates the sedimentation of solid particles in a viscous
fluid under gravity in two dimensions.

Keywords
--------
FEM, DEM, Generation

Description
-----------
The test demonstrates:
- How to generate a rectangular 2D mesh with physical boundary tags.
- How to define a gravity-driven sedimentation of dense particles.
- How to couple the particle solver (`scontact`) and the fluid solver in MigFlow.
- How to export derived fields (pressure, velocity, porosity) for post-processing.

.. code-block:: python
  
  import os, sys, shutil
  import numpy as np
  import gmsh
  from migflow import fluid, scontact, time_integration
  

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)
  

Geometrical parameters and mesh generation
------------------------------------------
A rectangular 2D mesh is generated using Gmsh with named boundaries.

.. code-block:: python
  
  height = 0.6  # domain height [m]
  width = 0.4  # domain width [m]
  mesh_size = 0.01  # mesh size [m]
  h = 1.5e-1  # particle bed height [m]
  w = 4e-1  # particle bed width [m]
  origin = [-width / 2, -height / 2]  # lower-left corner of the domain
  eps = 1e-8  # numerical tolerance
  
  
  # Function to generate the mesh
  def gen_mesh(width, height, mesh_size, origin=np.array([0, 0])):
      """Generate a rectangular 2D mesh with named boundaries."""
      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 [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(lambda dim, tag, x, y, z, lc: mesh_size)
      gmsh.model.mesh.generate(2)
  
  
  gen_mesh(width, height, mesh_size, origin)
  

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

.. code-block:: python
  
  g = np.array([0, -9.81])  # gravity [m/s²]
  r = 1.5e-3  # particle radius [m]
  rhop = 1500  # particle density [kg/m³]
  rho = 1000  # fluid density [kg/m³]
  nu = 1e-6  # kinematic viscosity [m²/s]
  mu = nu * rho  # dynamic viscosity [Pa·s]
  

Particle Problem Initialization
-------------------------------
Particles are placed in a rectangular patch at the top of the domain.

.. code-block:: python
  
  p = scontact.ParticleProblem(2)
  p.set_fixed_contact_geometry(0)
  p.load_msh_boundaries(None, ["Top", "Left", "Right", "Bottom"])
  

Grid Generation
---------------
Define particle positions in a grid with small random perturbations.

.. code-block:: python
  
  lx = w - 2 * r - eps
  ly = h
  y0 = height / 2 - r - eps
  step = 2 * r + eps
  x = np.arange(-lx / 2, lx / 2 - step, step) + step / 2
  n = x.shape[0]
  x_off = lx - n * step
  x += x_off / 2
  y = np.arange(y0, y0 - h, -step)
  x, y = np.meshgrid(x, y)
  for xi, yi in zip(x.flat, y.flat):
      p.add_particle((xi, yi), r, r**2 * np.pi * rhop)
  

Fluid Problem Initialization
----------------------------
The fluid solver is initialized using the generated mesh. All walls are no-slip.

.. code-block:: python
  
  f = fluid.FluidProblem(2, g, mu, rho)
  f.load_msh(None)
  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)
  
  
  def get_fields(fluid):
      """Return derived fields for output."""
      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),
      }
  
  

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

.. code-block:: python
  
  outf = 1  # number of iterations between outputs
  dt = 1e-2  # time step [s]
  tEnd = 2  # final time [s]
  t = 0
  i = 0
  

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

.. code-block:: python
  
  mass = np.pi * p.r() ** 2 * rhop
  while t < tEnd:
      print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
      # Coupling: update particle data in the fluid solver
      f.set_particles(
          p.delassus(),
          p.volume(),
          p.position(),
          p.velocity(),
          p.omega(),
          p.contact_forces(),
      )
      # Write outputs
      if i % outf == 0:
          f.write_mig(outputdir, t, get_fields(f))
          p.write_mig(outputdir, t)
  
      # Fluid and particle update
      f.implicit_euler(dt)
      fext = g * mass + f.compute_node_force()
      time_integration._advance_particles(p, fext, dt, 1, 1e-4 * r)
  
      # Time stepping
      t += dt
      i += 1
  

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

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

