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

Rotating Drum
=============
This example illustrates a circular shear flow generated by a rotating inner
cylinder (drum) surrounded by an outer stationary wall. Grains are placed in
the annular gap and interact with the fluid flow through drag and pressure force.

Keywords
--------
FEM, DEM, LMGC

Description
-----------
The test demonstrates:
- How to define moving boundary conditions as functions of position.
- How to generate an initial particle packing in a circular domain.
- How to couple the particle and fluid solvers in MigFlow.
- Optional use of `lmgc90` instead of `scontact` for contact resolution.

.. code-block:: python
  
  import os, sys, shutil, subprocess, time, random
  import numpy as np
  from migflow import fluid, scontact, time_integration
  
  use_lmgc90 = False
  if use_lmgc90:
      from migflow import lmgc90Interface
  

Output Directory and Mesh Generation
------------------------------------
The output directory is created, and a 2D annular mesh is generated using Gmsh.

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

Initial Particle Generation
---------------------------
Function to create an initial packing of grains within the annular domain.

.. code-block:: python
  
  def genInitialPosition(initdir, r, rout, rin, rhop, use_lmgc90):
      """Generate a particle packing and write it to an initial output file.
  
      Parameters
      ----------
      initdir : str
          Directory where the initial state is written.
      r : float
          Maximum particle radius.
      rout : float
          Outer radius of the drum.
      rin : float
          Inner radius of the drum.
      rhop : float
          Particle density.
      use_lmgc90 : bool
          If True, uses LMGC90 instead of scontact for contacts.
      """
      # Create particle problem
      p = scontact.ParticleProblem(2)
      print(mesh_filename)
      p.load_msh_boundaries(mesh_filename, ["Outer", "Inner"], material="Steel")
  
      # Define grid of possible particle centers
      x = np.arange(rout, -rout, 2.5 * -r)
      x, y = np.meshgrid(x, x)
      R2 = x**2 + y**2
  
      # Keep only points inside the annular region
      keep = (R2 < (rout - r) ** 2) & (R2 > (rin + r) ** 2)
      x = x[keep]
      y = y[keep]
  
      # Add particles with slight density variation
      for xi, yi in zip(x, y):
          if yi < rout:
              rhop1 = random.choice([rhop * 0.9, 1.1 * rhop, rhop])
              p.add_particle((xi, yi), r, r**2 * np.pi * rhop1, "Sand")
  
      p.write_mig(initdir, 0)
  
  

Physical and Numerical Parameters
---------------------------------

.. code-block:: python
  
  g = np.array([0, 0])  # no gravity
  rho = 1.253e3  # fluid density [kg/m³]
  rhop = 1000  # particle density [kg/m³]
  nu = 1e-3  # kinematic viscosity [m²/s]
  mu = rho * nu  # dynamic viscosity [Pa·s]
  rout = 0.0254  # outer radius [m]
  rin = 0.0064  # inner radius [m]
  r = 397e-6 / 2  # grain radius [m]
  

Particle Problem Initialization
-------------------------------

.. code-block:: python
  
  genInitialPosition(initdir, r, rout, rin, rhop, use_lmgc90)
  
  if use_lmgc90:
      friction = 0.1
      lmgc90Interface.scontactTolmgc90(initdir, 2, 0, friction)
      p = lmgc90Interface.ParticleProblem(2)
  else:
      p = scontact.ParticleProblem(2)
      p.read_mig(initdir, 0)
      p.set_friction_coefficient(0.1, "Sand", "Sand")
      p.set_friction_coefficient(0.1, "Sand", "Steel")
  p.set_fixed_contact_geometry(0)
  
  

Fluid Problem Initialization
----------------------------
The fluid solver is initialized, with a rotating inner boundary and fixed outer wall.

.. code-block:: python
  
  f = fluid.FluidProblem(2, g, nu * rho, rho)
  f.load_msh(mesh_filename)
  f.set_mean_pressure(0)
  f.set_wall_boundary("Outer", velocity=[0, 0])
  f.set_strong_boundary(
      "Inner",
      velocity=[lambda x: (x[:, 1] / rout) * 0.1, lambda x: (-x[:, 0] / rout) * 0.1],
  )
  
  # Coupling with particles
  f.set_particles(
      p.delassus(),
      p.volume(),
      p.position(),
      p.velocity(),
      p.omega(),
      p.contact_forces(),
  )
  

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

.. code-block:: python
  
  dt = 5e-3  # time step [s]
  tEnd = 10  # total simulation time [s]
  outf = 10  # number of iterations between outputs
  t = 0  # time [s]
  i = 0  # iteration counter

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

.. code-block:: python
  
  mass = np.pi * p.r() ** 2 * rhop  # particle masses
  while t < tEnd:
      print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
      if i % outf == 0:
          p.write_mig(outputdir, t)
          f.write_mig(outputdir, t)
      time_integration.iterate(
          f,
          p,
          dt,
          use_predictor_corrector=False,
          external_particles_forces=g * mass,
          contact_tol=1e-3 * r,
      )
      t += dt
      i += 1

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

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

